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CHAPTER  I 


INTRODUCTION 


Basic  queueing  theory  begins  with  an  arrival  process,  a 
service  mechanism,  and  a  queue  discipline.  Practical  appli¬ 
cations  can  extend  this  beginning  into  a  network  where  the 
nodes  would  be  service  mechanisms  of  one  or  more  servers. 
Such  applications  include  communication  systems,  computer 
time  sharing  processes,  medical  care  facilities,  assembly 
operations,  and  so  on.  The  analysis  of  such  networks 
involves  the  solution  of  large  scale  systems  of  equations 
and  computational  problems  of  large  dimensions.  Due  to  the 
intractability  of  the  mathematical  models,  computer  simula¬ 
tion  is  a  commonly  employed  analysis  approach. 

Simulation,  however,  is  an  experimental  approach  rather 
than  an  analytical  one,  and  presents  a  host  of  issues 
inherent  in  sampling.  These  issues  include  the  choice  of 
input  distributions,  statistical  methods  for  analyzing  out¬ 
put,  the  comparison  of  alternative  systems,  model  verifica¬ 
tion  and  validation,  and  techniques  used  to  improve  the 
precision  of  estimators.  The  last  issue  is  commonly 
referred  to  as  variance  reduction  and  is  the  topic  of  this 


research . 
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A  simulation  of  queueing  networks  is  partially  driven  by 
sampling  realizations  of  random  variables;  therefore,  the 
outputs  produced  are  also  random  variables.  These  outputs 
are  generally  mapped  into  estimates  of  interest  through  an 
output  function  (e.g.  a  sample  mean  for  instance).  These 
estimators  possess  sampling  distributions  usually  having 
unknown  means.  The  precision  of  these  estimators  is  meas¬ 
ured  by  their  variances:  the  smaller  the  variance  the 
greater  the  precision.  Therefore,  reducing  the  variance  is 
a  method  for  increasing  the  efficiency  of  the  simulation. 

One  technique  employed  for  variance  reduction  is  the  use 
of  control  variates.  This  technique  uses  the  correlation 
between  specified  random  variables  to  achieve  a  variance 
reduction.  One  type  of  control  variate  is  the  external 
control,  which  is  obtained  by  simulating  a  similar  system 
whose  performance  measures  can  be  analytically  computed  or 
closely  approximated.  A  variance  reduction  can  be  obtained 
if  the  output  of  the  second  simulation  is  positively  corre¬ 
lated  with  its  counterpart  from  the  original  simulation. 

A  number  of  queueing  networks  can  be  categorized  as 
Jackson  networks,  for  which  the  analytical  computation  of 
various  performance  measures  is  possible.  Jackson  networks 
have  been  considered  as  possible  control  variates  for  simu¬ 
lating  more  general  networks.  There  are  at  least  two  ways 
in  which  Jackson  networks  can  be  used  to  obtain  control 
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variates.  The  first  method,  as  described  above,  involves 
running  a  second  simulation  of  a  similar  Jackson  network 
and  using  the  corresponding  output  of  this  second  simula¬ 
tion  as  the  control  variate.  This  method  is  commonly 
referred  to  as  external  control  variates.  A  formidable 
drawback  with  this  approach  is  the  cost  of  the  second  simu¬ 
lation. 

The  other  method  for  obtaining  a  control  variate  is  to 
use  the  difference  between  two  performance  measures  calcu¬ 
lated  from  the  Jackson  model  as  the  control  variate.  One 
measure  is  computed  by  substituting  the  known  input  parame¬ 
ters  into  the  Jackson  equations .  These  parameters  could  be 
the  mean  arrival  or  service  rates  used  to  drive  the  simula¬ 
tion.  The  other  performance  measure  is  computed  by  substi¬ 
tuting  estimates  of  these  same  parameters  obtained  from  the 
simulation  into  the  Jackson  equations.  This  type  of  con¬ 
trol  variate  is  referred  to  as  an  analytic  control  since  it 
is  obtained  from  an  analytical  operation  rather  than  a  sec¬ 
ond  simulation.  The  advantage  of  this  approach  is  eliminat¬ 
ing  the  cos t  of  the  second  simulation.  This  is  a  new 
approach. 

The  purpose  of  this  research  is  to  study  the  effective¬ 
ness  of  Jackson  networks  as  external  and  analytic  controls 
for  queueing  network  simulations.  The  approach  taken  is  to 
experiment  with  a  small  but  representative  set  of  networks 
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with  an  eye  toward  drawing  general  conclusions  about  the 
performance  of  this  variance  reduction  technique.  Nelson 
[15,16]  notes  that  prior  knowledge  of  the  system  in  ques¬ 
tion  is  a  key  component  in  the  selection  of  an  appropriate 
variance  reduction  technique.  The  conclusions  drawn  from 
this  research  should  provide  the  analyst  some  prior  knowl¬ 
edge  for  selecting  the  appropriate  variance  reduction  tech¬ 
nique. 

This  research  will  attempt  to  add  to  this  prior  knowl¬ 
edge  by  studying  the  performance  of  Jackson  based  external 
and  analytic  controls  on  various  queueing  performance  meas¬ 
ures,  by  investigating  the  impact  of  the  service  distribu¬ 
tions,  traffic  intenstity,  and  network  structure.  In 
addition,  the  suitability  of  automating  this  approach  and 
areas  of  future  research  will  be  discussed.  The  remainder 
of  this  work  includes  a  background  on  queueing  networks  and 
control  variates,  the  methodology  used  in  this  research, 
results  and  conclusions. 


CHAPTER  II 


BACKGROUND 


The  purpose  of  this  chapter  is  to  present  an  integrated 
review  of  the  literature  relevant  to  this  research.  The 
review  is  divided  into  three  sections:  the  first  presents  a 
brief  introduction  to  queueing  networks  and  formally 
defines  a  Jackson  network;  the  second  presents  the  theory 
and  development  of  control  variates,  and  the  third  discuss¬ 
es  the  results  of  the  control  variate  techniques  applied  to 
queueing  network  simulations. 

QUEUEING  NETWORKS 

In  general,  queueing  networks  are  classified  as  open  or 
closed  networks.  In  an  open  network  customers  arrive  from 
outside  the  network;  this  characteristic  is  called  exoge¬ 
nous  arrivals.  In  general  customers  may  enter  the  network 
at  any  node.  Customers  then  proceed  through  the  network 
according  to  their  needs  or  in  some  random  manner  and  may 
depart  the  network  from  any  node.  A  closed  network  is  simi¬ 
lar  in  structure  to  an  open  network;  however,  there  are  no 
exogenous  arrivals  and  customers  never  depart  the  network. 


There  is  always  some  fixed  number  of  customers  present  in  a 


closed  network.  Figure  1  shows  examples  of  the  open  and 
closed  network  types .  This  network  contains  a  number  of 
points  where  customer  routing  decisions  must  be  made.  These 
points  are  called  switches  and  their  operation  is  governed 
by  switch  rules.  These  rules  may  be  imposed  externally  to 
the  system  (e.g.  a  routing  or  dispatch  form),  internal  to 
the  system  (e.g.  the  server  at  node  1  may  determine  whether 
a  customer  goes  to  node  2  or  3),  or  the  rules  may  be  deter¬ 
mined  by  the  customers  (e.g.  customer  selects  the  shortest 
waiting  line). 


External 

Arrival 


-> 

Departure 


Closed  Network (Number  of  Customers*  C) 


Figure  1 :  Open  and  Closed  Queueing  Networks 
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While  the  systems  in  Figure  1  illustrate  the  idea  of 
open  and  closed  networks,  more  detailed  symbols  are  needed 
to  model  more  complicated  structures.  Consider  the  network 
in  Figure  2. 


Feedback 


External 

Arrival 


Decomposition 

Switch 


Recomposition  Departure 
Switch 


Sample  Queueing  Network 


There  are  two  basic  types  of  switches:  decomposition  and 
recomposition.  A  decompostion  switch  splits  a  single  stream 
of  customers  into  a  number  of  streams.  A  recomposition 
switch  merges  a  number  of  streams  into  one  superposed 


stream. 

Another  possible  feature  is  feedback.  A  feedback  point 
is  one  where  customers  may  be  directed  to  repeat  a  service 
node;  direction  is  provided  through  feedback  rules. 

There  are  three  principal  methods  for  analyzing  queueing 
networks:  first,  analyze  the  network  as  a  whole;  second, 
decompose  the  network  into  subnetworks;  and  third. 


use  computer  simulation.  In  his  survey  paper  Disney  £33 
notes  that  the  difficulties  in  mathematically  analyzing 
queueing  networks  arise  from  flow  properties  rather  than 
physical  properties.  Once  customers  enter  a  network  the 
combinatorial  effect  of  service  mechanisms,  switch  and 
feedback  rules,  and  queue  disciplines  alter  the  flow  within 
the  system  for  the  individual  customer. 

Many  of  the  techniques  for  analyzing  networks  as  a  whole 
are  based  on  the  research  of  J.R.  Jackson  [7,83*  The  major 
thrust  in  this  area  has  been  studying  the  queue  length  pro¬ 
cess  and  most  of  the  known  results  are  for  steady-state 
behavior.  The  primary  obstacles  encountered  are  finding  the 
solutions  of  large  scale  systems  of  equations. 

The  second  approach,  decomposition,  attempts  to  break 
the  network  down  into  subnetworks  whose  characteristics  are 
well  known.  The  most  commonly  used  point  for  decomposition 
is  at  the  switches.  There  are  two  basic  technical  problems 
with  this  approach:  first,  determining  the  effect  of 
switching  rules  on  the  stochastic  properties  of  network 
flow;  and  second,  determining  the  result  of  recombination 
of  the  subnetworks.  The  primary  obstacles  encountered  are 
probabilistic  as  opposed  to  algebraic,  and  involve  computa¬ 
tional  problems  of  large  dimension. 

The  third  approach,  computer  simulation  is  probably  the 
most  commonly  employed  for  general  networks.  Simulation  is 
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an  experimental  approach  rather  than  an  analytical 
approach.  The  obstacle  faced  is  how  to  analyze  the  measures 
obtained  from  the  simulation.  Often  simulation  output  is 
used  to  estimate  a  population  mean.  In  general,  the  output 
is  correlated  and  highly  variable.  Estimation  and  the  con¬ 
trol  of  the  variance  of  estimators  is  important,  and  is 
reflected  in  the  validity  and  width,  respectively,  of 
interval  estimators  of  these  popluation  means. 

The  characteristics  of  a  queueing  network  are  principal¬ 
ly  determined  by  the  arrival  processes,  service  mechanisms, 
queue  disciplines,  switches,  and  feedback  rules.  The  model 
formulated  by  Jackson  [73  properly  defines  these  character¬ 
istics  so  as  to  facilitate  a  generalization  of  the  M/M/s 
queue  (Kendall  notation  meaning  exponential  interarrival 
times/  exponential  service  times/s  servers)  to  an  intercon¬ 
nected  open  network  of  service  nodes.  The  defining  charac¬ 
teristics  of  a  Jackson  network  are  listed  below: 

1.  The  network  contains  more  than  one  service  node. 

2.  Each  node  can  be  a  single  or  multiple  server  queue 
with  each  server  having  an  identical  exponential  ser¬ 
vice  time  distribution.  Service  times  are  independent. 

3.  Arrivals  from  outside  the  network  occur  in  a  Poisson 
fashion.  Outside  arrivals  to  any  node  are  independent. 

4.  Arrivals  at  any  given  node  may  come  from  outside  the 
network  or  from  any  other  nodes. 


5. 


The  effective  arrival  rate  at  every  node  is  less  than 
its  potential  service  rate. 

6.  When  a  customer  completes  service  at  a  node,  he  may 
leave  the  network  or  be  routed  to  another  node. 

7.  There  is  unlimited  waiting  space  at  every  service 
node . 

8.  The  queue  discipline  is  first  come  first  served. 
Although  the  parameters  are  fairly  restrictive,  the  mod¬ 
el  is  still  quite  general.  Subsets  of  Jackson  networks  are 

1.  A  finite  number  of  M/M/s  queues  in  tandem;  tandem  net¬ 
works  have  only  one  exogenous  arrival  point  and  one 
path  through  the  network. 

2.  An  acyclic  network  of  M/M/s  queues;  these  are  networks 
where  customers  may  visit  a  node  only  once. 

3.  A  network  of  M/M/s  queues  with  feedback. 

Jackson  proved  the  important  result  that  in  steady 
state  conditions,  each  node  in  the  network  functions  as  an 
independent  M/M/s  queue  with  Poisson  input.  This  fact 
allows  the  decomposition  into  subnetworks  and  the  pursuant 
application  of  M/M/s  results.  Another  important  fact  is 
that  although  feedback  destroys  the  Poisson  property  of  the 
input  stream,  the  nodes  of  the  Jackson  network  continue  to 
function  as  though  the  input  process  was  Poisson. 

Lemoine  [13]  presented  a  survey  of  equilibrium  results 
for  general  Jackson  networks.  If  the  network  is  open,  the 
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equilibrium  rate  of  flow  through  node  i,  e^  ,  is  the  sum  of 
the  external  input  rate,  A^,  and  the  total  rate  of  inter¬ 
nal  inputs  to  node  i.  This  balance  equation  can  be  written 
as 

e±«  A±  +  £  e j  i=l,...,N  (1) 

where  r^  is  the  probability  a  customer  is  routed  from 
node  j  to  node  i,  and  N  is  the  number  of  service  nodes  in 
the  network. 

Since  the  effective  arrival  rate  at  each  node  i  must  be 
less  than  its  potential  service  rate,  or  else  the  number  of 
customers  in  the  system  will  continually  grow  as  time  goes 
on,  the  traffic  intensity  must  satisfy 


9i=ei/(siui)<l  i=l . N  (2) 

where  s^  is  the  number  of  servers  and  u^  is  the  service 
rate  at  node  i . 

Another  equilibrium  flow  condition  derived  from  the  open 
network  is  that  the  total  input  flow  rate  must  equal  the 
external  departure  rate.  For  any  node  x,  the  probability 
that  any  customer  leaves  the  network  is 


(3) 
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therefore. 


m 


IK- 1 1  * 


1=1  *  i=i 


In  his  work  Jackson  [7]  used  as  a  state  variable  a  vec¬ 
tor  whose  components  represent  the  number  of  customers 
present  at  each  node  in  the  network.  His  analysis  showed 
that  under  equilibrium  conditions  the  probabilty  of  the 
system  being  in  state  k,  p(k),  could  be  factored  into  a 
product  of  the  marginal  probabilities: 

P(k^ ,  •  • .  , kN ) ~P^  1^^)  ...  Pjj ( k^ )  ( 5 ) 
where 


si'1 


P±(0)  = 


(e./Ui)k  +  (ei/ui)  1 

k-T  s±l  (1-^) 


p.(0)  (ei/ui)J 
ki 


k— 0 ,1, • • • , 


Pi (k)=< 


Pi (0 )  ( e±/u± ) 


k>Si 


I 

so  (s.) 

l  '  i' 


m 

S 

% 
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The  preceding  two  equations  can  be  recognized  as  those 
of  the  basic  M/M/s  model  with  effective  arrival  rate,  e^, 
replacing  A  This  generalization  permits  the  decomposi¬ 
tion  of  the  Jackson  network  into  a  collection  of  multi¬ 
server  subnetworks. 

Using  the  above  results  and  Little's  [14]  formula  the 
following  measures  of  performance  can  be  obtained: 


i 

Long  run  queue  length  (LQi)=pi(0)  (ei/ui)  (pA) 


si  (1-  ); 


Long  run  node  length  (L^Jslq^  +  e^/u^ 


i 


The  long  run  node  length  is  the  sum  of  the  number  of  cus¬ 


tomers  in  service  and  the  number  in  the  queue. 
Long  run  queue  time  (WQ^)=LQi/e^ 


(10) 


Long  run  node  time  (W.)*WQ  +  1/u. 

i  l  '  i 


(11) 


Nelson  [17]  extended  the  results  for  Jackson  networks 
by  deriving  the  probability  distribution  for  the  total 
waiting  time  (excluding  service  time)  for  a  customer  to 
pass  completely  through  the  network.  This  result  was 
obtained  by  the  convolution  of  waiting  time  distributions 
at  each  node. 
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The  sojourn  time  of  a  customer  in  a  network  is  the  time 
spent  at  ea\:h  of  the  nodes  visited  (queue  and  service  time 
combined)  plus  travel  time  between  the  nodes.  Travel  time 
in  a  Jackson  network  is  assumed  to  be  zero.  For  most  net¬ 
works  the  sojourn  time  problem  is  unsolved.  Burke  [1]  and 
Reich  [18]  present  results  for  small  special  case  networks. 

Gordon  and  Newell  [6]  analyze  a  closed  network  of  N 
interconnected  nodes  and  C  customers.  Each  node  has  s.  . 

l  ' 

i=l,...,N,  parallel  servers  each  with  service  rate  u^.  The 
routing  from  node  to  node  is  the  same  as  a  Jackson  networx 
except  customers  do  not  depart  the  network.  The  system  is 
equivalent  to  some  open  networks  where  the  number  of  cus¬ 
tomers  cannot  exceed  C.  The  authors*  principal  result  was 
an  expression  for  the  equilibrium  distribution  at  each 
node.  The  expression  is  factored  into  product  terms  for 
each  node  with  the  exception  of  an  unknown  normalizing  con¬ 
stant  that  reflects  the  interaction  between  nodes. 

Buzen  [2]  developed  an  iterative  technique  for  deter¬ 
mining  the  normalizing  constant.  He  also  derived  the  margi¬ 
nal  distributions  of  the  number  of  customers  present  at  the 
nodes,  the  expected  number  of  customers  at  each  node  and 
the  steady  state  utilizations. 

Solberg  [19]  developed  a  computationally  efficient  meth¬ 
od  for  computing  the  normalizing  constant. 
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In  summary,  the  analysis  of  Jackson  networks  have  the 
following  limitations: 

1.  Service  time  distributions  must  be  exponential. 

2.  Service  nodes  must  have  identical  servers. 

3.  Only  probabilistic  routing  between  nodes  is  permitted. 

4.  Customer  oriented  performance  measures  such  as  sojourn 
times  are  difficult  to  obtain  for  other  than  special 
cases . 

5.  Travel  time  between  nodes  is  assumed  to  be  zero. 

The  scarcity  of  analytical  results  for  other  than  special 
cases  makes  using  network  models  difficult  for  practical 
applications .  In  general  computer  simulation  often  becomes 
the  analysis  approach  and,  as  mentioned  previously,  estima¬ 
tors  of  the  performance  measures  will  possess  some  degree 
of  variability.  Reducing  this  variability  to  increase  the 
estimator's  precision  becomes  a  major  concern.  One  way  of 
addressing  this  concern  is  the  use  of  control  variates. 


CONTROL  VARIATES 

The  central  idea  of  control  variates  is  to  use  the  correla¬ 
tion  between  specified  random  variables  to  achieve  a  vari¬ 
ance  reduction.  A  random  variable,  C,  is  a  control  variable 
for  the  random  variable  Y,  if  it  has  a  known  expectation, 
*Y,  and  is  correlated  with  Y. 


Let  Y  be  an  unbiased  estimator  of  ©,  the  quantity  of 
interest,  obtained  from  a  single  simulation  run.  Then  for 


1  . 


any  constant  b,  an  estimator  of  Y  can  be  written  as 


Y(b)  =  Y  -  b(C  -  J  ) 


(12) 


Equation  (12)  is  also  an  unbiased  estimator  of  0.  The  vari¬ 


ance  of  Y(b)  is  given  by 


Var[Y(b)]=Var[Y]  +  b^Var[C]  -  2bCov[Y,C] 


which  is  the  same  as 


Var[Y(b) ]=Var[Y]  +  b2Var[C]  -  2b?Var[Y]  Var[c] 


(14) 


where  is  the  coefficient  of  correlation  between  Y  and  C. 


The  value  of  b,  b  ,  which  minimizes  the  Var[Y(b)J  can  be 


found  by  differentiating  (13)  with  respect  to  b  and  is  giv¬ 


en  by 


b  =  Cov[Y,C] 


(15) 


Var[C] 


Substituting  the  above  into  (12)  yields  the  optimal  control 


variate  estimator  Y(b  ).  The  variance  of  this  estimator  is 


then 


Var[ Y ( b  )]=VarCY]  (1-9^) 


Equation  (16)  indicates  the  greater  the  correlation  between 
Y  and  C,  the  greater  the  reduction  in  variance. 

Kleijnen  [9]  discusses  extensions  to  multiple  control 
variates 

Y(b)=  Y  -  £  U?> 

where  n  is  the  number  of  control  variables. 

Law  and  Kelton  [12]  present  two  general  methods  for 
obtaining  control  variables.  The  first  is  to  use  input  ran¬ 
dom  variables,  such  as  arrival  rates,  service  rates,  and 
routing  probabilities,  since  their  expectations  are  known 
and  the  sign  of  the  correlation  with  the  output  may  be 
known.  This  type  of  control  variate  is  known  as  .internal  or 
concomitant.  Since  they  are  generated  by  the  simulation  to 
obtain  the  outputs,  using  them  adds  little  to  the  cost  of 
the  simulation. 

A  second  method  for  obtaining  control  variates  is  to 
simulate  a  similar  system  whose  desired  performance  measure 
can  be  analytically  computed.  This  simulation  uses  the  same 
random  numbers  as  the  first  simulation  to  induce  positive 
correlation.  The  corresponding  output  of  the  second  simu¬ 
lation  can  then  be  used  as  the  control  variate.  This  type 
of  control  variate  is  called  an  external  control  variate. 
The  desired  outcome  is  that  the  output  of  the  second  simu¬ 
lation  is  positively  correlated  with  its  counterpart  from 


the  original  simulation.  Unlike  internal  control  variates 
the  cost  of  a  second  simulation  is  incurred,  which  in  some 
cases  may  be  prohibitive.  Thus  the  covariance  between  the 
outputs  will  have  to  be  larger  than  for  the  internal  con¬ 
trol  variates  to  make  this  approach  worthwhile. 

A  third  method  for  obtaining  control  variates,  suggested 
by  Nelson  £15],  is  an  amalgam  of  the  internal  and  external 
approaches.  He  suggests  simulating  the  system  to  obtain  the 
desired  performance  measures  and  the  means  of  the  input 
parameters  observed  during  the  simulation  run.  The  control 
variate  is  derived  by  substituting  these  observed  input 
means  into  a  parametric  analytical  model  of  a  similar  sys¬ 
tem.  The  mean  of  this  control  variate  would  be  derived  in  a 
similar  fashion,  except  the  known  input  means,  rather  than 
the  observed  means,  would  be  substituted  into  the  paramet¬ 
ric  model.  Expressed  in  the  linear  control  variate  format, 
the  control  estimator  of  Y  would  be 


Y(b)=  Y  -  b  (A 


-  r 


(18) 


where  Y  is  the  crude  estimator  obtained  from  the  simula- 

—  4 

tion;  A=g(X^),  is  the  control  variate  where  is  the 
observed  mean  of  the  input  X^^  driving  the  simulation  that 
produced  Y  and  i=l,...c,  where  c  is  the  number  of  input 
parameters  in  the  parametric  model  function  g;  and 
^"=g(E[X<])  where  EfX.3  is  the  input  mean.  Equation  (18) 
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need  not  be  unbiased  since  the  expectation  of  a  function  is 
not  in  general  a  function  of  expectations. 

From  a  cost  standpoint  the  analytic  approach  has  an 
advantage  over  the  true  external  in  that  the  cost  of  the 
second  simulation  is  avoided.  The  effectiveness  of  this 
approach  using  the  Jackson  network  as  the  parametric  model 
is  the  focal  point  of  this  research. 

Once  a  control  variate  method  has  been  selected  the 
problem  of  specifying  the  control  coefficient,  b,  must  be 
addressed. 

Consider  the  case  where  there  is  only  one  control  vari¬ 
ate,  C,  and  (12)  is  used  as  the  control  estimator.  Then  the 
optimal  value  of  b,  b*,  is  expressed  in  (15).  In  general 
Cov[Y,C]  and  the  Var[C}  are  not  known;  therefore,  b*  needs 
to  be  estimated. 

Kleijnen  [9]  presents  a  method  for  estimating  b*  from 
the  simulation  results.  He  suggests  replacing  Cov[Y,Cj  and 
VarEeJ  with  their  sample  equivalents.  Consider  making  n 
independent  replications  to  obtain  n  independent  and  iden¬ 
tically  distributed  (iid)  observations  of  Y  and  C.  Then  the 
sample  covariance  of  Y  and  C,  Cov[Y,C],  is  given  by 


n 


CovEY,  C>_J_  Y  (Y.-Y)(C,-C) 
n-1  1  1 


(19) 


i=l 


and  the  sample  variance  of  C,  V&r[Cj  is 


n 


V&r[C>_L  ]T  (cA-c)2 

n-1 


(20) 


i=l 
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then  the  estimator  for  b,  b  is  given  by 


b  =  Cov[Y,C] 


(21) 


V&r[C] 


This  produces  a  final  point  estimator  of  © 


Y(b* )=  Y-b*(C-7  ) 


(22) 


It  should  be  noted  that  Y($>  )  may  not  be  unbiased  since  $>* 
and  C  are  not  usually  independent,  since  £>*  is  a  function 
of  C  as  given  by  (21).  The  author  discusses  two  tech¬ 
niques,  splitting  and  jackknifing,  for  reducing  the  bias  of 
Y(b* ) . 

The  case  of  multiple  control  variates  is  addressed  by 
Lavenberg  and  Welch  [10j.  The  following  notation  is  adopted 
to  rewrite  (17)  in  matrix  form.  Let  JL  be  a  column  vector, 
andJL'  be  its  transpose.  Then  is  a  column  vector  of  Q 
control  variates  and  X  is  the  mean  vector  corresponding  to 
C  where  Let  b  be  a  vector  of  constants.  Then  an 
estimator  of  ©  is 


'(£)=  Y  -  h  (£  -  X  ) 


The  vector  b  which  minimizes  Var[Y(ji)]  it 
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where  Y.  c  is  the  covariance  matrix  of  C  and  i£yC  is  the 

Q-dimensional  vector  whose  components  are  the  covariances 

between  Y  and  the  C^,  i=l,...,Q.  This  leads  to  a  minimum 

* 

variance  for  Y(Ji  ): 

Var[Y(fc*)>  U“RyC)  Var[Y]  (25) 

where 


Rvc  -  —  YC  ^C1  —VC  (26> 

Var[Y] 

2  2 
and  (1-Ryc)  is  called  the  minimum  variance  ratio.  RyC 

is  the  square  of  the  multiple  correlation  coefficient 

between  Y  and  C. 

* 

As  with  the  single  control  variate  b  is  unknown  and 
must  be  estimated.  An  estimator  of  b*  is 


*  * 


I-i  <r 

C  YC 


(27) 


where  is  the  sample  covariance  matrix  and  _££yC  is  the 
sample  covariance  vector. 

To  derive  interval  estimates  the  authors  consider  obser¬ 
vations  from  J  statisically  independent  but  otherwise  iden¬ 
tical  runs.  Then  £  would  be  a  vector  of  control  variates 
whose  components  are  the  values  of  £^  on  the  jth  replica¬ 


tion.  Then 


22 


Y.-k*(C.-2)  (28) 

and 

J 

?(b*)=_L  y  Y.(b*)  (29) 

j  b,  3 


_  A  * 

In  general  Y(b  )  is  not  an  unbiased  estimator  of  6  and  the 
t-distribution  with  J-l  degrees  of  freedom  cannot  be  used 
to  derive  the  interval  estimate.  The  authors  derive  con¬ 
fidence  intervals  for  the  multiple  control  case  based  on 
the  assumption  that  the  vector  (Y.C^, . . . ,cQ)  has  a  multi¬ 
variate  normal  distribution.  Under  this  assumption  stan¬ 
dard  regression  techniques  can  be  used  to  produce 
var[Y(b*)] 


The  above  equation  indicates  that  if  J,  the  number  of  rep¬ 
lications  ,  is  not  large  with  respect  to  Q,  the  number  of 
control  variates,  the  variance  reduction  produced  by 
(l-RyC)  will  be  diminished.  The  authors  report  experi¬ 
mentation  which  showed  this  factor  accurately  predicted 
losses  in  variance  reduction. 

APPLICATION  OF  CONTROL  VARIATES 

The  control  variate  approach  was  applied  to  queueing  net¬ 
work  simulations  by  Lavenberg,  Moeller  and  Welch  [11],  and 
Gaver  and  Schedler  [53*  A  summary  of  these  works  follows. 

Lavenberg  et  al.  £11 3  considered  the  application  of 
internal  control  variates  to  a  broad  class  of  closed  net¬ 
works.  These  networks  allowed  priorities,  blocking,  differ¬ 
ent  customer  types  and  arbitrary  service  time  distribu¬ 
tions.  Their  network  consisted  of  n  finite  interconnected 
nodes  with  one  or  more  servers,  and  d=l,...,D  customer 
types.  A  type  d  event  is  the  departure  of  a  type  d  custom¬ 
er.  Customer  routing  through  the  network  was  controlled  by 
an  (nxn)  transition  probability  matrix.  The  following  meas¬ 
ures  were  obtained:  W^(d),  the  expected  queue  time  for  type 
d  customers  at  node  i;  A(d),  the  expected  rates  at  which 
events  occur  for  type  d  customers;  T(d),  the  expected  time 
for  type  d  customers  to  cycle  through  the  network  and 
return  to  the  first  node. 


The  authors  experimented  with  three  types  of  control 
variates,  all  of  which  were  internal  control  variates.  The 
first,  work  variables,  represented  the  sum  of  service  times 
for  type  d  customers  at  node  i  for  type  d  events  in  the 
system.  The  second,  flow  variables,  represented  the  frac¬ 
tion  of  type  d  events  at  node  i.  The  third,  service  vari¬ 
ables,  represented  the  sample  mean  service  times  for  type  d 
customers  at  node  i. 

The  authors  reported  substantially  larger  variance 
reductions  using  work  variables  as  opposed  to  flow  or  ser¬ 
vice  variables.  Experimentation  was  then  limited  to  work 
variables.  For  a  network  consisting  of  four  to  six  nodes 
and  one  customer  type,  they  report  predicted  actual  vari¬ 
ance  ratios  using  six  control  variates  (Q=6)  of  .30  to  .85. 
Predicted  actual  variance  ratios  were  obtained  by  multiply¬ 
ing  the  estimated  minimum  variance  ratios  by  the  theoreti¬ 
cal  loss  factor.  Estimated  variance  ratios  ranged  from  .16 
to  .77,  and  are  ratios  of  the  variance  of  the  point  estima¬ 
tor  with  work  variables  to  the  variance  of  the  crude  point 
estimator.  The  largest  variance  reductions  for  waiting 
times  were  achieved  at  the  node  having  the  largest  utiliza¬ 
tion  factor. 

Wilson  and  Pritsker  [21]  performed  a  similar  study  using 
standardized  work  variables.  Work  variables  are  standard¬ 
ized  for  a  specific  time  period  by  correcting  each  variable 


observed  by  its  mean  and  standard  deviation.  This  is  per¬ 
formed  so  the  control  variates  would  be  asymptotically  sta¬ 
ble,  ensuring  efficiency  gains  are  sustained  over  increas¬ 
ing  statistic  accumulation  intervals.  The  authors  report 
variance  reductions  of  20  to  90  percent.  They  stated  their 
standardized  work  variables  could  not  be  extended  to  simu¬ 
lations  of  open  and  mixed  networks. 

Gaver  and  Schedler  £5]  applied  external  control  vari¬ 
ates  to  a  closed  network.  Their  study  was  the  only  one 
found  reporting  results  for  external  controls.  Their  net¬ 
work  contained  two  service  nodes  each  offering  three  dif¬ 
ferent  types  of  services.  They  allowed  for  priority  service 
and  a  mixture  of  arbitrary  and  exponential  service  time 
distributions.  Steady  state  utilization  factors  were  the 
performance  measures  of  interest.  Their  control  variate 
was  the  utilization  obtained  from  the  simulation  of  a  simi¬ 
lar  but  numerically  tractable  model. 

Results  were  reported  for  control  variate  estimators 
using  a  control  coefficient  equal  to  one  and  an  estimated 
optimal  control  coefficient  based  on  (21).  For  the  control 
coefficient  equal  to  one  variance  reductions  of  51  to  99 
percent  were  achieved  with  one  exception:  a  node  with  99 
perecent  utilization  produced  a  29  percent  increase  in  var¬ 
iance.  For  the  estimated  optimal  control  coefficient  case 
variance  reductions  of  81  to  99  percent  were  achieved.  The 
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authors  note  the  latter  estimates  may  not  be  unbiased  since 
b  was  estimated  from  the  data;  however,  this  bias  decreases 
as  the  sample  size  increases.  No  direct  conclusions  could 
be  drawn  about  the  relationship  of  utilization  and  variance 
reduction.  The  results  did  indicate  a  trend  in  which  uti¬ 
lization  estimates  with  large  sample  variances  showed  the 
largest  variance  reductions  after  the  application  of  con¬ 
trol  variates. 

As  stated  in  Chapter  1  the  purpose  of  this  research  is 
to  study  the  effectiveness  of  Jackson  networks  as  external 
and  analytic  control  variates  for  open  queueing  networks. 
Figure  3  places  this  research  in  the  context  of  previous 
work  in  this  area. 


Control  Variates  in  Queueing  Network  Simulations 


Closed 


I- Internal  Lavenberg  et  al 


Open(this  research) 

•External 


+: 


*■  External  Gaver  and  Schedler  LAnalytic 
Figure  3:  Context  of  the  Research 

The  study  of  external  and  analytical  control  variates 
applied  to  open  queueing  networks  is  largely  without  prec¬ 
edent.  The  network  structures  to  be  tested  and  the 


methodology  of  carrying  out  these  tests  are  very  experimen¬ 
tal.  This  makes  it  difficult  to  predict  in  advance  the 
suitability  of  these  controls  for  this  class  of  simulation 
problems . 


The  primary  objective  of  this  research  is  to  investigate 
the  effectiveness  of  the  Jackson  model  as  a  control  variate 


for  queueing  network  simulation.  Three  different  network 
structures  were  investigated,  each  meeting  the  restrictions 
of  the  Jackson  model  with  the  exception  of  service  time 
distributions.  Service  distributions  investigated  were  the 
exponential,  the  Weibull,  and  the  uniform.  For  each  network 
two  types  of  control  variates  were  obtained:  the  tradition¬ 
al  external  control  variate  and  the  analytic  control  vari¬ 
ate.  These  controls  were  obtained  to  estimate  the  steady 
state  measures  of  server  utilization  factor  and  customer 
queue  time  at  each  node. 

The  utilization  factor  was  selected  because  it  serves  as 
an  indicator  of  the  level  of  activity  or  degree  of  con¬ 
gestion  at  a  particular  node  i.  The  queue  time  provides  the 
long  run  waiting  time  a  customer  will  experience  in  a  given 
queue  (excluding  service  time),  and  when  applied  in  Lit¬ 
tle's  formula,  WQ^,  yields  the  long  run  queue 
length.  Additionally,  the  steady  state  values  for  the 


total  time  spent  at  a  service  node  i,  W^,  and  the  number  of 
customers  at  the  node,  ,  can  be  found  by  substituting  WQ^ 
and  LC^  into  (9)  and  (11). 

TYPES  OF  CONTROL  VARIATES 

The  classical  external  control  variate  requires  two  simula¬ 
tions.  The  first  simulation  is  of  the  network  of  interest 
and  estimates  the  utilization  factor  and  queue  time  for 
each  node.  The  second  simulation  is  of  a  Jackson  network 
approximation  of  the  original  system.  Since  the  exponential 
distribution  yields  the  Jackson  model  itself,  external  con¬ 
trol  variates  are  obtained  only  for  the  Weibull  and  uniform 
cases.  For  each  of  these  distributions  a  second  simulation 
was  run  with  common  random  numbers  using  the  means  of  the 
Weibull  and  uniform  distributions  as  the  parameters  of  the 
exponential,  and  the  two  desired  performance  measure  esti¬ 
mates  were  obtained.  Using  these  means  and  the  Jackson  mod¬ 
el  equations,  the  corresponding  steady  state  measures  were 
obtained  analytically.  The  control  variate  estimators  at 
each  node  i  (i=l,...,N)  for  the  utilization  factor,  RO.C^, 
and  queue  time,  WQ.C^,  based  on  external  control  variates 
are  given  by 
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where  RO.S^  is  the  crude  estimator  of  the  utilization  fac¬ 
tor  at  node  i  obtained  from  the  simulation;  RO.E^  is  the 
external  control  variate  obtained  from  the  second  simula¬ 
tion;  and  RO(J)^  is  the  analytic  value  of  the  steady  state 
Jackson  network  based  on  the  parameters  A  ,  u^,  _r  ;  WQ.S^ 
WQ.E^,  and  WQ(J)^  are  defined  similarly  for  the  queue 
times . 

The  major  drawback  of  this  type  of  control  variate  is 
the  cost  of  the  second  simulation  and  the  associated  prob¬ 
lem  of  synchronization.  A  system  is  synchronized  when  a 
random  number  used  for  a  purpose,  such  as  arrival  or  ser¬ 
vice  times,  in  one  system  is  used  for  the  same  purpose  in 
\ 

the  other  systems  being  compared.  The  random  numbers  are 
those  generated  from  the  uniform  (0,1)  distribution  and 
mapped  through  an  inverse  transformation  into  the  desired 
probability  distribution,  such  as  Poisson  or  exponential. 
Synchronization  tries  to  solve  the  problem  of  insuring  that 
differences  between  the  two  simulations  are  due  to  model 
performance  and  not  random  number  sequences  or  coding 
structure.  If  the  systems  were  not  synchronized  the  com¬ 
parison  of  control  variate  performance  might  be  influenced 
by  the  misapplication  of  random  numbers. 

In  contrast,  the  analytic  control  variate  requires  only 
simulation  of  the  network  of  interest.  To  obtain  the  ana¬ 
lytic  control  variate  additional  coding  is  added  to  the 
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simulation  program  to  record  the  vectors  of  observed  mean 
arrival  rates,  A.  •  service  rates,  u.  ,and  the  observed  frac¬ 
tion  of  customers  routed  to  the  various  nodes,  r  .The  ana¬ 
lytic  control  variates  at  each  node  i  for  the  utilization 
factor,  RO ( J )  ^ ,  and  queue  time,  WQ(J)i,  are  obtained  from 
the  Jackson  model  equations 


RO(J)i=ei/(siui)  (34) 

WQ(J)i=LQi/ei  (35) 

where  e^  is  defined  in  (1)  and  LQ^  is  defined  in  (8).  The 
analytic  Jackson  values  for  the  steady  state  utilization 
factors  and  queue  times  were  calculated  in  the  same  way  as 
the  external  control  variates.  Analytic  control  variate 
estimates  could  then  be  calculated  from 


RO.Ci=RO.Si-b(RO(J)i-RO(J)i)  (36) 

WQ.Ci*WQ.Si-b(WQ(J)i-WQ(J)i)  (37) 

Equations  (32)  and  (33)  are  the  same  as  (36)  and  (37)  with 
the  exception  of  the  control  variate  terms.  In  (32)  and 
(33)  the  control  variates  RO.E^,  and  WQ.E^^  are  obtained 
from  a  second  simulation  of  a  network  modelled  as  a  Jackson 
network.  In  (36)  and  (37)  the  control  variates  are 
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A  yy  ^ 

obtained  by  substituting  .A,  u,  r  into  the  Jackson  equa¬ 


tions  . 


The  tradeoff  in  using  the  analytic  control  variate  is 


A  A  A 

the  additional  code  required  to  obtain  A «  u,  _r  .  This 


additional  coding  is  insignificant  relative  to  the  cost  of 


a  second  simulation.  Once  the  controls  or  data  needed  to 


compute  the  controls  has  been  obtained  the  computational 


effort  to  obtain  the  control  variate  estimates  is  the  same 


for  both  type  of  control  variates. 


One  drawback  of  the  analytic  control  variate  is  at  high 


traffic  intensities  (  $--9)  it  is  possible  to  obtain  sample 
values  for  ju ,  £  which  violate  Jackson  model  assump¬ 


tions,  specifically  ei/(siui)<l.  This  does  not  permit  the 


calculation  of  an  analytic  control  variate  based  on  the 


Jackson  equations.  A  possible  solution  could  be  to  observe 


the  effective  arrival  rates,  e.  ,  rather  than  observing 


and  .£.  This  approach  was  employed  with  one  of  the  net¬ 


works.  Further  studies  of  these  two  calculation  methods  is 


required  to  determine  the  benefits  and  tradeoffs  of  each 


method.  Another  drawback  is  that,  in  general,  the  expecta¬ 


tion  of  a  function  is  not  equal  to  the  function  of  the 


expectation  (e.g.  E[R0(J)^3^H0(J)i) ;  however,  because  it  is 


a  method  of  moments  estimator,  it  is  consistent.  Therefore 


using  the  observed  mean  arrival  and  service  rates  in  the 


Jackson  model  equations  will  result  in  a  biased  control 
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variate.  A  study  to  determine  the  severity  of  this  bias  is 
an  area  open  to  future  reseach,  since  reduced  variance  at 
the  expense  of  significantly  large  mean  squared  error  is 
unacceptable. 


NETWORKS 

To  obtain  a  representative  appraisal  of  the  effectiveness 
of  Jackson  model  control  variates,  three  networks  with  dif¬ 
ferent  structure  and  complexity  were  simulated  using  common 
random  numbers.  The  first  network  consists  of  two  nodes, 
each  with  its  own  external  arrival  process.  Customers  com¬ 
pleting  service  at  each  node  may  be  routed  to  the  other 
node  for  service  or  may  depart  the  system  entirely. 

The  second  network  is  a  three  node  tandem,  acyclic  net¬ 
work.  Tandem  means  the  nodes  are  arranged  in  series  and 
acyclic  means  customers  will  visit  each  node  once.  External 
arrivals  occur  only  at  the  first  node  where  customers  com¬ 
plete  service  and  move  to  the  second  and  then  third  nodes 
for  service.  Departure  from  the  network  occurs  only  when 
service  is  completed  at  the  third  node. 

The  third  network  consists  of  four  nodes  with  an  exter¬ 
nal  arrival  process  at  the  first  node.  Customers  completing 
service  at  the  first  node  are  routed  either  to  the  second 
or  third  nodes  and  then  on  to  the  fourth.  Customers  com¬ 
pleting  service  at  the  fourth  node  may  be  fed  back  to  the 
first  node  or  depart  the  system  entirely. 


Three  service  time  distributions  were  studied  for  each 
of  the  above  networks.  The  first  distribution,  the  exponen¬ 
tial,  is  the  requisite  for  the  Jackson  model.  It  has  an 
infinite  tail,  is  highly  variable,  and  provides  analytical¬ 
ly  tractable  performance  measures  for  comparing  the  control 
variate  estimators.  The  second  distribution,  the  Weibull, 
is  similar  to  the  exponential.  It  also  has  an  infinite 
tail,  but  it  is  not  as  variable  as  the  exponential,  as 
characterized  by  the  coefficient  of  variation.  By  setting 
the  Weibull  shape  parameter  to  2a  humped  distribution  was 
obtained,  thereby  providing  another  reference  point  to 
measure  the  effectiveness  of  the  Jackson  controls.  The 
third  distribution,  the  uniform,  was  selected  for  its  mark¬ 
ed  difference  from  the  exponential.  It  has  finite  range  and 
is  considerably  less  variable.  The  selection  of  these  dis¬ 
tributions  provides  three  references  for  studying  the  Jack- 
son  controls:  the  exponential,  the  requisite  for  the  Jack- 
son  model,  infinite  in  the  tail,  and  highly  variable;  the 
Weibull,  similar  to  the  exponential  but  humped  in  our  exam¬ 
ples;  and  the  uniform,  a  finite  range  distribution  with 
considerably  less  variability. 

Other  features  of  the  networks  studied  include  traffic 
intensity  and  the  number  of  servers  at  each  node.  To  study 
the  effect  of  congestion  on  control  variate  performance 
both  high,  Q  =.90,  and  low,  Q=.50,  traffic  intensities 
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were  studied.  Traffic  intensity  measures  the  fraction  of 
the  systems  service  capacity  being  utilized  on  the  average 
by  arriving  customers.  Traffic  intensities  close  to  1  mean 
there  will  rarely  be  idle  servers,  so  customers  will  be 
found  backing  up  into  the  queues.  Queue  times  and  lengths 
will  therefore  be  larger.  Single  and  multiple  servers  were 
studied  in  each  network. 

In  summary,  the  basic  experiment  was  to  investigate  two 
types  of  Jackson  control  variates,  external  and  analytic, 
for  estimating  the  utilization  factors  and  waiting  times  in 
three  different  queueing  networks.  The  control  variates  for 
each  network  were  obtained  for  three  service  time  distribu¬ 
tions:  exponential,  Weibull,  and  uniform;  and  at  both  high 
and  low  traffic  intensities.  Figure  4  provides  an  outline 
of  the  basic  experiment  for  a  given  network.  Figures  5-7 
provide  schematics  and  parameters  for  each  of  the  three 
networks . 
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ANALYTIC  ANALYTIC  EXTERNAL  ANALYTIC  EXTERNAL 

Figure  4:  Outline  of  the  Basic  Experiment 
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Figure  5:  Network  I  and  parameters 
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Figure  6:  Network  II  and  parameters 
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Figure  7:  Network  III  and  parameters 


EXPERIMENTAL  DESIGN 

The  basic  computational  steps  required  to  obtain  the  ana¬ 
lytic  control  variate  estimators  are  defined  by  (36)  and 
(37). 

Previously  it  was  mentioned  that  under  high  traffic 
intensity  conditions  it  was  not  always  possible  to  obtain  a 
control  variate  for  a  given  batch  (defined  below).  In  these 
cases  the  vectors  X. ,  _u_,  _r  may  produce  utilization  factors 
greater  than  or  equal  to  one,  a  violation  of  the  Jackson 
model  assumptions.  To  handle  this  case  the  ratio  of  arrival 
rate  to  service  rate  was  set  equal  to  .9999  whenever  the 
utilization  factor  was  greater  than  or  equal  to  one.  This 
in  effect  is  the  use  of  control  variates  from  all  the 
batches  at  the  expense  of  introducing  some  bias  into  the 
control  estimator. 

The  discussion  pertaining  to  (21)  outlined  the  proce¬ 
dure  for  obtaining  the  estimated  optimal  control  coeffi¬ 
cient.  The  procedure  required  n  independent  replications  to 
obtain  n  iid  observations  of  the  crude  estimator  and  its 
control.  From  these  n  values  the  sample  variance  and  covar¬ 
iance  terms  of  (21)  can  be  calculated. 

Since  the  performance  measures  of  interest  are  steady 
state  measures,  an  initial  bias  period  for  each  replication 
would  have  to  be  deleted.  If  the  initial  transient  is  long, 
as  will  be  discussed  later,  the  cost  of  deleting  n  initial 


bias  periods  becomes  excessive.  To  avoid  this  costly 
approach  a  simulation  run  consisting  of  J  approximately 
independent  batches  of  time  length  t  was  used  instead  of 
the  n  independent  replications. 

Following  this  procedure,  the  Y  and  C  of  equations 
(19)  through  (21)  are  now  replaced  by  Ko.S^,  and  1*0(3)^, 
where  Eo.S^  is  the  batch  mean  for  the  crude  estimator  of 
the  utilization  factor  at  node  i  in  batch  j,  where 
i=l,...,N  and  j=l,...,J.  R0(3)^_.  is  the  analytic  control 
variate  derived  from  u.  »  JL  in  batch  j.  Let 


RO.S .  .=  _1 
1  J 


i=l 


i  N 


(38) 


KO(J)  •  .  =  _1 

^  T 


I 


j  =  l 


KD(J)  .  . 
13 


(39) 


Var[RO(J) . ]=  1 
1  J-l 


J 

]T  (RO(J)ij-RO(j)i.  )2 

j=l 


(40) 


(41) 


Cov[RO.S . ,RQ( J) . ]=  1 
1  1  J-l 


J 

]T  (RO.Sij-RO.Si.  ) 
j  =  1 

x(RO(J)i;.-RO(J)i.  ) 


Using  (21)  the  estimated  value  for  the  optimal  control 
coefficient  of  the  utilization  factor  at  node  i  is 


b*==  Cov[RO.Si<RO(J)i]  (42) 

An  analagous  procedure  was  followed  to  obtain  /b*  for  WQ.C. 

*  * 

The  value  of  b  computed  in  (42)  can  be  used  to  compute 
the  control  variate  estimate  for  the  run  by  using 

RO.Ci=RO.Si.-b*(RO(J)i.-RO(J)i)  (43) 

Since  the  utilization  factor  is  a  time  persistent  per¬ 
formance  measure,  computing  RO.S^  from  batches  of  equal 
time  length  produces  an  unbiased  estimator  assuming  each  of 
the  are  identically  distributed.  This  is  not  the 

case  for  the  queue  time  performance  measure.  WQ  is  a 


Var[RO(( 


discrete  performance  measure,  therefore  batching  by  time 
produces  a  random  number  of  customer  queue  times  observed 
in  each  of  the  j  batches.  The  overall  mean  queue  time  for  a 


run,  WQ.S^.  is  given  by 


(44) 


times  during 

D 

WQ.S.  .*  J.  Y  WQ.S.  (45) 

i  •jj  m 

n=l 

where  D  is  the  total  number  of  queue  times  observed  in  the 
run.  To  accommodate  the  discrete  case,  the  grand  mean  for 
all  the  queue  times  for  the  run,  wQ.S^,  was  used  to  calcu¬ 
late  the  control  variate  estimate.  W&.S^  is  given  by 


WQ 


U 

•si-4  I  TO-Sii 


j=l 

J 


I 


WQ.Sin 


j  ncB  . 
J  D 


where  B_.  is  the  set  of  all  indices  of  queue 
((j-l)t,jt),  and  d_.=  |Bj(  .  Therefore 


The  control  variate  estimator  of  the  queue  time  for  the 
run  is  given  by 

WQ.Ci=WQ.Si-b*{WQ(J)i.-WQ(J)i)  (47) 

The  identical  approach  is  taken  to  obtain  the  external 
control  variate  estimators.  The  only  difference  being  sub¬ 
stituting  KO.E^  and  WQ.E^  for  RO(J)^  and  WQ(J)^  respective¬ 
ly  in  (43)  and  (47).  In  practice  these  external  controls 
would  be  obtained  by  simulating  the  network  of  interest  as 
a  Jackson  network.  Since  exponential  service  yields  the 
Jackson  model  itself,  the  values  of  RO.E^  and  WQ.  E^  equal 
RO. and  WQ.S i  from  the  networks  with  exponential  service 
times . 

The  batch  means  approach  serves  two  purposes.  First, 

the  J  batches  per  run  provide  a  sequence  of  observations  to 

a  * 

compute  b  and  the  control  variate  estimators.  Second,  K 
runs  of  J  batches  each  can  be  obtained  by  simulating  a 
total  of  K-J  batches.  This  will  produce  a  sequence  of  K 
control  variate  estimates  so  that  the  properties  of  the 
estimator  can  be  evaluated.  The  primary  design  issues  with 
this  approach  are  the  batch  length  t,  and  the  number  of 
batches  to  be  collected  within  a  run,  and  within  the  entire 


experiment. 


The  number  of  batches  selected  for  a  particular  run  was 
based  on  cost  considerations  and  the  loss  of  variance 

A  * 

reduction  caused  by  estimating  b  .  This  loss  was  expressed 
in  (31)  as  a  function  of  the  number  of  control  variates  and 

A  * 

the  number  of  batches  used  to  estimate  b  . 

While  multiple  controls  are  possible,  this  research 
studies  only  a  single  commensurate  control;  that  is,  the 
corresponding  performance  measure  for  the  Jackson  network. 
The  single  control  approach  was  adopted  to  contain  the  cost 
of  gathering  control  variate  statistics  and  to  facilitate 
automation.  In  addition,  if  the  number  of  control  variates 
is  large  with  respect  to  the  number  of  batches,  considera¬ 
ble  loss  in  variance  reduction  will  result.  Since  Q=1  the 
loss  factor,  LF,  can  be  expressed  as 

LF  =  (48) 

J-3 

where  J  is  the  number  of  batches  in  a  run.  Table  1  lists 
various  loss  factors  and  their  corresponding  number  of 
batches . 

Based  on  the  above  comparison  and  cost  factors,  the  number 

of  batches,  J,  was  set  at  25.  The  tradeoff  of  estimating 
*  * 

b  was  then  a  5  percent  loss  in  variance  reduction. 

The  batch  length  issue  centers  on  choosing  a  time 
length,  t,  large  enough  to  secure  approximate  independence 
between  the  batch  means.  It  is  assumed  the  output  sequence 
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Table  1 


Loss  Factor  Comparison 


J  Loss  Factors (LF) 


10 

1.14 

15 

1.08 

20 

1.06 

25 

1.05 

30 

1.04 

50 

1.02 

of  crude  estimators  and  their  counterpart  control  variates 
is  covariance  stationary,  and  the  batch  length  will  be 
large  enough  so  that  the  resulting  batch  means  will  be 
approximately  normally  distributed.  To  select  the  batch 
length,  t,  an  independence  test  given  by  Fishman  [4]  was 
employed  for  each  network  to  evaluate  the  independence  of 
queue  times  at  each  node. 

The  results  of  this  testing  produced  a  batch  size  and 
corresponding  number  of  batches  based  on  a  type  I  error 
level  of  .05.  The  results  for  the  node  in  each  network 
requiring  the  largest  number  of  observations  per  batch  are 
listed  in  Table  2. 

The  number  of  batches  was  fixed  at  25  for  cost  and  loss 
factor  considerations  as  previously  discussed.  Since  the 


Table  2 


Results  for  Batch  Means  Independence  Test 


NETWORK 

NODE 

NO.  BATCHES 

SIZE 

TOTAL  OBS. 

1 

1 

49 

64 

3136 

2 

1 

23 

128 

2944 

3 

4 

12 

512 

6144 

This  assumption  facilitates  the  conversion  from  discrete 
batch  size  to  continuous  time  batch  length.  This  is  done  by 
dividing  the  total  number  of  observations  from  Table  2  by 
25.  Results  are  listed  in  Table  3. 


Table  3 


Selected  Batch  Lengths 


Selected 


NETWORK 

OBS. /Batch 

Batch  Length (time 

units ) 

1 

125.4 

150 

2 

117.8 

200 

3 

245.8 

300 

The 

major  concern 

in 

deciding  the  number  of 

runs  K,  the 

"macro” 

replications 

for 

computing  the  point 

and  interval 

estiamtes,  was  cost.  The  CPU  time  required  to  simulate  the 
two  node  network  for  10  runs  was  1.3  seconds,  the  time  for 
20  runs  was  1.7  seconds,  an  approximate  24  perecnt  increase 


in  CPU  time.  This  time  will  also  increase  with  network 


research  and  its  associated  cost,  the  number  of  runs,  K, 
for  each  experiment  was  fixed  at  10. 

The  performance  measures  being  investigated  are  steady 
state  means;  therefore,  a  procedure  to  eliminate  the  ini¬ 
tial  transient  was  employed  at  the  start  of  each  experi¬ 
ment.  To  approximate  the  length  of  the  initial  transient, 
a  pilot  run  listing  the  cumulative  mean  of  the  queue  time 
at  each  node  in  intervals  of  100  time  units  was  executed. 
The  results  showed  that  from  empty  and  idle  conditions  the 
build  up  to  steady  state  was  very  slow.  The  system  appeared 
very  erratic  during  the  first  10,000  time  units  before  set¬ 
tling  down  in  a  more  predicable  region  around  the  steady 
state  conditions.  Therefore  a  conservative  policy  of  elim¬ 
inating  the  statistics  collected  during  the  first  10,000 
time  units  after  starting  from  empty  and  idle  conditions 
was  adopted. 

The  pilot  runs  indicated  that  those  networks  with  high 
feedback  tended  to  reach  steady  state  sooner  than  those 
with  lesser  or  no  feedback.  We  can  only  speculate  that  the 
high  feedback  tends  to  congest  the  system  sooner,  which  in 
turn  has  a  stabilizing  effect  allowing  the  system  to  reach 
steady  state  at  a  faster  rate.  It  should  be  noted  that 
this  effect  was  studied  only  at  high  traffic  intensities 
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The  results  of  the  control  variate  experimentation  on  each 
of  the  three  networks  are  listed  in  this  chapter.  Three 
control  variate  estimates  are  reported:  the  analytic,  the 
modified  analytic,  and  the  external.  The  experimentation 
was  conducted  to  produce  control  variate  estimates  using 
control  coefficients  equal  to  1  and  equal  to  the  estimated 
optimal  control  coefficient  $>*.  When  the  control  coeffi¬ 
cient  was  set  to  1  the  variance  reductions  for  the  utiliza¬ 
tion  factor  estimates  were  slightly  greater  than  those 

A  * 

obtained  for  b  ;  however,  this  was  not  true  for  the  queue 
time  estimates.  Variances  of  these  estimates  were  greatly 
increased  when  the  control  coefficient  was  set  to  1;  vari- 

A  it 

ance  reductions  for  this  measure  were  achieved  only  when  b 

was  used.  Therefore  results  are  reported  only  for  esti- 
*  * 

mates  based  on  b  . 

The  effectiveness  of  a  particular  control  variate  will 
be  reported  in  terms  of  the  variance  reduction  ratio;  that 
is,  the  ratio  of  the  variances  of  the  control  variate  esti¬ 
mate  to  the  crude  estimate.  Ratios  greater  than  1  indicate 


5 


an  increase  in  the  estimate's  variance.  A  90  percent  confi¬ 
dence  interval  is  computed  for  each  ratio.  Ratios  involv¬ 
ing  the  expenditure  of  computer  effort  are  also  possible, 
but  not  considered  here. 

The  chapter  is  divided  into  three  sections  corresponding 
to  the  three  networks  studied.  Results  for  each  network  are 
reported  in  the  following  format:  a  table  listing  the 
steady  state  Jackson  values  for  each  performance  measure,  a 
table  listing  the  crude  estimator  and  its  variance  for  a 
given  traffic  intensity,  and  a  table  for  each  of  the  three 
control  variate  estimators  listing  the  point  estimate,  its 
variance,  the  variance  reduction  ratio,  and  the  upper  and 
lower  bounds  of  the  90  percent  confidence  interval.  Vari¬ 
ance  reduction  ratios  appear  under  the  ratio  column,  the 
lower  bound  under  the  L  column,  and  the  upper  bound  under 
the  U  column.  Estimates  with  variances  less  than  .00005  are 
reported  as  "<.00005". 

The  simulation  was  coded  using  SIMSCRIPT  II. 5  and  all 
required  output  written  to  a  file.  A  FORTRAN  program  was 
used  to  perform  the  control  variate  analysis.  Samples  of 
both  programs  are  in  the  Appendix. 


RESULTS  FOR  NETWORK  I 

Table  4  lists  the  steady  state  Jackson  values  for  the  uti¬ 
lization  factors  and  queue  times  of  network  I  depicted  in 
Figure  5. 


Table  4 


Steady  State  Jackson 

Values  for 

Network 

9  =-9 

5 

NODE 

RO  { J ) 

WQ(  J ) 

RO  ( J ) 

WQ(J) 

1 

.  9000 

7.4771 

.5000 

.4616 

2 

.8999 

9.1999 

.5000 

.4001 

The  following  table  lists  the  crude  estimates  and  their 
variances  for  each  of  the  three  distributions  tested  at  the 
high  traffic  intensity  level. 


Table  5 

Crude  Estimates  for  Network  I  (P-.9) 
Exponential  Service 


NODE 

RO.S 

VAR 

WQ.S 

VAR 

1 

.8872 

.0008 

5.9196 

2.1802 

2 

.9076 

.0006 

9.6711 

4.3639 

Weibull 

Service 

1 

.8934 

.0005 

4.1025 

.8210 

2 

.9033 

.0003 

4.9167 

1.0489 

Uniform 

Service 

1 

.8991 

.0002 

3.7309 

.5230 

2 

.9023 

.0002 

3.9522 

.4220 

Table 

6  lists 

results  for 

the  analytic 

control 

estimates 

at  the 

high  traffic 

intensity. 

Table  6  listed  results  for  analytic  control  variates 

^  A  A 

based  on  A.,  u,  and _£.  Additional  experimentation  was  con- 

a  A 

ducted  to  determine  if  other  combinations  of  _u,  and  r_ 
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Table  6 


Analytic  Estimates  for  Network  I  (P=.9) 
Exponential  Service 


NODE 

RO.C 

VAR 

L 

RATIO 

U 

1 

.9025 

.0001 

.0487 

.1550 

.4929 

2 

.9130 

.0002 

.1019 

.3240 

1.0303 

NODE 

WQ.C 

VAR 

L 

RATIO 

U 

1 

5.0694 

1.1656 

.1682 

.  5350 

1.7013 

2 

8.4098 

3.3497 

.2420 

.7695 

2.4470 

Weibull  Service 

NODE 

RO.C 

VAR 

L 

RATIO 

U 

1 

.9010 

.0001 

.0401 

.1275 

.4055 

2 

.9034 

<.00005 

.0471 

.1499 

.4767 

NODE 

WQ.C 

VAR 

L 

RATIO 

U 

1 

3.7813 

.5469 

.2095 

.6661 

2.1182 

2 

4.1912 

.5425 

.  1626 

.5172 

1.6447 

Uniform  Service 

NODE 

RO.C 

VAR 

L 

RATIO 

U 

1 

.9023 

<.00005 

.0718 

.2284 

.7263 

2 

.9017 

< .00005 

.0752 

.2392 

.7607 

NODE 

WQ.C 

VAR 

L 

RATIO 

U 

1 

3.3113 

.3808 

.2290 

.7281 

2.3154 

2 

3.4828 

.  3020 

.2250 

.7154 

2.2750 

would  produce  a  better  control  variate.  A  pilot  run  of  Net¬ 
work  I  showed  a  modified  analytic  control  variate  based  on 

A 

A.*  iL,  and  £•  that  is  based  on  the  observed  mean  arrival 
rates,  the  input  mean  service  time,  and  the  observed  rout¬ 
ing  probabilities,  was  the  most  promising.  These  modified 
analytic  control  variate  estimators  for  the  utilization 
factor  and  queue  time  are  denoted  as  RO(M)  and  WQ(M) 
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respectively.  Results  for  this  modified  estimator  are 
reported  in  Table  7. 


Table  7 


Modified  Estimates  for  Network  I  (p  =  .9) 
Exponential  Service 


NODE 

RO  (M) 

VAR 

L 

RATIO 

U 

1 

.8889 

.0005 

1974 

.6278 

1.9964 

2 

.9016 

.0004 

2145 

.6821 

2.1691 

NODE 

WQ(M) 

VAR 

L 

RATIO 

U 

1 

5.6553  1 

.6770 

2419 

.7692 

2.4461 

2 

9.3293  4 

.8245 

3477 

1.1055 

3.5155 

Weibull 

Service 

NODE 

RO  (M ) 

VAR 

L 

RATIO 

U 

1 

.8936 

.0002 

1398 

.4445 

1.4135 

2 

.8979 

.0002 

1504 

.4783 

1.5210 

NODE 

WQ(M) 

VAR 

L 

RATIO 

U 

1 

3.8313 

.5795 

2219 

.7058 

2.2444 

2 

4.3954 

.  5804 

1740 

.5533 

1.7595 

Uniform 

i  Service 

NODE 

RO  (M ) 

VAR 

L 

RATIO 

U 

1 

.8990  < 

.00005 

0716 

.2277 

.7241 

2 

.8970 

.0001 

0856 

.2722 

.8656 

NODE 

WQ(M) 

VAR 

L 

RATIO 

U 

1 

3.3723 

.3727 

2241 

.7126 

2.2661 

2 

3.4059 

.2790 

2079 

.6611 

2.1023 

Table  8  lists  results  for 

the 

external 

control  in  Net- 

work  I 

at  a  traffic 

intensity 

of  . 

9. 

Tables  9-12  list 

the  same 

type 

results 

for  Network  I  at 

a  traffic  intensity 

of  .5. 
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Table  8 


External 

Estimates 

for  Network  I 

(P  =  .9) 

Weibull 

Service 

NODE 

RO.C 

VAR 

L  RATIO 

U 

1 

.9017 

<.00005 

0146  .0464 

.1476 

2 

.8990 

.0001 

1397  .4443 

1.4129 

NODE 

WQ.C 

VAR 

L  RATIO 

U 

1 

4.6927 

.5279 

2022  .6430 

2.0447 

2 

4.6992 

.4128 

1237  .3935 

1.2513 

Uniform 

i  Service 

NODE 

RO.C 

VAR 

L  RATIO 

U 

1 

.9016 

.0001 

2048  .6514 

2.0715 

2 

.9004 

.0001 

1608  .5114 

1.6263 

NODE 

WQ.C 

VAR 

L  RATIO 

U 

1 

3.9298 

.8453 

5083  1.6164 

5.1402 

2 

3.8757 

.2829 

2108  .6704 

2.1319 

Table  9 

Crude  Estimates  for  Network  I  (P 

=  .5) 

Exponential  Service 

NODE 

RO.S 

VAR 

WQ.S 

VAR 

1 

.4940 

.0003 

.4481 

.0015 

2 

.  5038 

.0002 

.3922 

.0042 

Weibull 

Service 

1 

.4965 

.0002 

.2933 

.0003 

2 

.  5026 

.0001 

.2512 

.  0004 

Uniform  Service 

1 

.4986 

.0001 

.2420 

.0002 

2 

.5012 

.0001 

.2064 

.0002 

-.-w 


Table  10 


Analytic  Estimates  for  Network  I  (P=.5) 


Exponential  Service 


NODE 

RO.C 

VAR 

L 

RATIO 

U 

1 

.4996 

<.00005 

.0003 

.0009 

.0029 

2 

.4982 

<.00005 

.0011 

.0036 

.0114 

NODE 

WQ.C 

VAR 

L 

RATIO 

U 

1 

.4513 

.0005 

.  1058 

.3363 

1.0694 

2 

.3574 

.0018 

.1318 

.4190 

1.3324 

Weibull  Service 

NODE 

RO.C 

VAR 

L 

RATIO 

U 

1 

.4996 

<.00005 

.0003 

.0008 

.0025 

2 

.4982 

<.00005 

.0018 

.0056 

.0178 

NODE 

WQ.C 

VAR 

L 

RATIO 

U 

1 

.2950 

.0002 

.2254 

.7169 

2.2797 

2 

.2380 

.0003 

.2083 

.6625 

2.1068 

Uniform  Service 

NODE 

RO.C 

VAR 

L 

RATIO 

U 

1 

.4995 

< .00005 

.0003 

.0011 

.0035 

2 

.4981 

< . 00005 

.0024 

.0076 

.0242 

NODE 

WQ.C 

VAR 

L 

RATIO 

U 

1 

.2404 

.0002 

.2308 

.7341 

2.3344 

2 

.2017 

.0002 

.2558 

.8134 

2.5866 

Table  11 


Modified  Estimates  for  Network  I  (P=.5) 
Exponential  Service 


RO  (M ) 

VAR 

L 

RATIO 

U 

4913 

.0001 

.1152 

.3663 

1.1648 

4955 

.0001 

.1473 

.4685 

1.4898 

WQ(M) 

VAR 

L 

RATIO 

U 

4343 

.0007 

.  1419 

.4512 

1.4348 

3733 

.0031 

.2279 

.7246 

2.3042 

Weibull  Service 

RO  (M ) 

VAR 

L 

RATIO 

U 

4940 

<.00005 

.0796 

.2530 

.8045 

4940 

< .00005 

.0924 

.2939 

.9346 

WQ(M) 

VAR 

L 

RATIO 

U 

2854 

.0001 

.0796 

.2530 

.8045 

2357 

.0003 

.1998 

.6353 

2.0203 

Uniform  Service 

RO  (M ) 

VAR 

L 

RATIO 

U 

4970 

<.00005 

.0215 

.0685 

.2178 

4942 

< . 00005 

.0106 

.0338 

.1075 

WQ(M) 

VAR 

L 

RATIO 

U 

2362 

.0002 

.2685 

.8538 

2.7151 

1958 

.0002 

.2869 

.9124 

2.9014 

Weibull  Service 


NODE 

RO.C 

VAR 

L 

RATIO 

U 

1 

.5007  < 

.  00005 

.0392 

.  1246 

.  3962 

2 

. 5007  < 

.00005 

.0379 

.1206 

.3835 

NODE 

WQ.C 

VAR 

L 

RATIO 

U 

1 

.3025 

.  0004 

.3276 

1.0418 

3.3129 

2 

.2552 

.0005 

.3473 

1.1044 

3.5120 

Uniform  Service 

NODE 

RO.C 

VAR 

L 

RATIO 

U 

1 

.5008  < 

.  00005 

.2564 

.8155 

2.5933 

2 

.4998  < 

.00005 

.1793 

.  5703 

1.8136 

NODE 

WQ.C 

VAR 

L 

RATIO 

U 

1 

.2492 

.0006 

.8307 

2.6416 

8.4003 

2 

.2061 

.0002 

.2659 

.8455 

2.6887 

RESULTS 

FOR  NETWORK 

II 

Table  13 

lists  the 

steady 

state 

Jackson 

results 

utilization  factors  and  queue  times  of  Network  II  depicted 
in  Figure  6. 


Table  13 

Steady  State  Jackson  Values  for  Network  II. 


?  =-9 

II 

NODE 

RO  ( J ) 

WQ(  J ) 

RO  ( J ) 

WQ(  J ) 

1 

.9001 

8.1089 

.  5000 

.  5000 

2 

.8999 

7.6667 

.  5000 

.  3333 

3 

.  9001 

7.3466 

.  5000 

.2368 
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Tables  14-17  list  results  for  Network  II  at  a  traffic 


intensity  of  .9. 


Table  14 


Crude  Estimates  for  Network  II  (P  = . 9 ) 


Exponential  Service 


NODE 

1 

2 

3 


RO.S 

.8872 

.8947 

.8906 


VAR 

.0005 

.0001 

.0004 


WQ.S 

6.3363 

7.3719 

6.1086 


VAR 

1.1857 

4.4346 

3.0807 


Weibull  Service 


.8941 

.8939 

.8918 


.0002 

.0001 

.0002 


4.2966 

2.9552 

2.2526 


.4055 

.4286 

.2985 


Uniform  Service 


.8964 

.8946 

.8950 


,0001 

.0001 

,0001 


4.1722 
1.0114 
.  5184 


.3479 

.0388 

.0072 


Tables  18-21  list  results  for  Network  II  at  a  traffic 


intensity  of  . 5 . 


•NV.YsV 


MB 


y  v 


Table  15 


Analytic  Estimates  for  Network  II 


Exponential  Service 


NODE 

RO.C 

VAR 

L 

RATIO 

U 

1 

.9008 

<.00005 

.0239 

.0759 

.2414 

2 

.9010 

< . 00005 

.  1134 

.3607 

1.1470 

3 

.8985 

.0001 

.1013 

.3222 

1.0246 

NODE 

WQ.C 

VAR 

L 

RATIO 

U 

1 

5.6212 

.2973 

.0788 

.2507 

.7972 

2 

6.4935 

1.8028 

.1278 

.4065 

1.2927 

3 

5.6702 

2.5672 

.2608 

.8295 

2.6378 

Weibull  Service 

NODE 

RO.C 

VAR 

L 

RATIO 

U 

1 

.8997 

< . 00005 

.0197 

.0625 

.1988 

2 

.8999 

<.00005 

.0398 

.1267 

.4029 

3 

.8978 

< . 00005 

.0619 

.1968 

.6258 

NODE 

WQ.C 

VAR 

L 

RATIO 

U 

1 

3.8640 

.1924 

.1492 

.4744 

1.5086 

2 

2.7459 

.3320 

.2436 

.7747 

2.4635 

3 

2.1985 

.2648 

.2790 

.8872 

2.8213 

Uniform  Service 

NODE 

RO.C 

VAR 

L 

RATIO 

U 

1 

.9015 

< . 00005 

.0477 

.1516 

.4821 

2 

.8995 

<.00005 

.0371 

.1181 

.3756 

3 

.8996 

< . 00005 

.0481 

.1530 

.4865 

NODE 

WQ.C 

VAR 

L 

RATIO 

U 

1 

3.8384 

.4355 

.1775 

.5645 

1.7951 

2 

.9515 

.0314 

.2549 

.8105 

2.5774 

3 

.4910 

.0054 

.2373 

.7545 

2.3993 

:A 


•  »*■  rw  v*  v  \r 


H 


fr 

r* 
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Table  16 


Modified 

Estimates 

for  Network  II 

(P  =  .9) 

Exponential  Service 

NODE 

RO  (M ) 

VAR 

L 

RATIO 

U 

1 

.8805 

.0003 

.2214 

.  7041 

2.2390 

2 

.8971 

.0001 

.3334 

1.0602 

3.3714 

3 

.8917 

.0004 

.2730 

.8682 

2.7609 

NODE 

WQ(M) 

VAR 

L 

RATIO 

U 

1 

6.1244 

1.1639 

.3087 

.9816 

3.1215 

2 

7.1038 

3.7865 

.2685 

.8539 

2.7154 

3 

6.0812 

3.4835 

.3539 

1.1255 

3.5791 

Weibull  Service 

NODE 

RO  (M ) 

VAR 

L 

RATIO 

U 

1 

.8932 

.0001 

.1515 

.4818 

1.5321 

2 

.8961 

< . 00005 

.2073 

.6592 

2.0963 

3 

.8934 

.0001 

.1563 

.4971 

1.5808 

NODE 

WQ(M) 

VAR 

L 

RATIO 

U 

1 

4.0569 

.3096 

.2401 

.7634 

2.4276 

2 

2.8120 

.3552 

.2607 

.8289 

2.6359 

3 

2.2232 

.2692 

.2836 

.9020 

2.8684 

Uniform  Service 

NODE 

RO  (M ) 

VAR 

L 

RATIO 

U 

1 

.8988 

<.00005 

.1184 

.3766 

1.1976 

2 

.8965 

< . 00005 

.0599 

.1906 

.6061 

3 

.8970 

<.00005 

.0476 

.1513 

.4811 

NODE 

WQ(M) 

VAR 

L 

RATIO 

U 

1 

3.9215 

.2612 

.2361 

.7507 

2.3872 

2 

.9486 

.0314 

.2549 

.8105 

2.5774 

3 

.4953 

.0056 

.2435 

.7743 

2.4623 

Table  17 

External  Estimates  for  Network  II  ( 


& 


:w 


Weibull  Service 

RO.C 

VAR 

L 

RATIO 

U 

.9009 

< • 00005 

.0195 

.0619 

.1968 

.8976 

<.00005 

.  1520 

.4835 

1.5375 

.8972 

< . 00005 

.0265 

.0843 

.2681 

WQ.C 

VAR 

L 

RATIO 

U 

5.2104 

.4523 

.3508 

1.1154 

3.5470 

3.1890 

.1155 

.0847 

.2694 

.8567 

2.5559 

.1728 

.1821 

.5790 

1.8412 

Uniform  Service 

RO.C 

VAR 

L 

RATIO 

U 

.9003 

.0001 

.2906 

.9242 

2.9390 

.8965 

.0001 

.3332 

1.0596 

3.3695 

.8968 

< . 00005 

.1527 

.4855 

1.5439 

WQ.C 

VAR 

L 

RATIO 

U 

4.4618 

.3898 

.3523 

1.1203 

3.5626 

1.0130 

.0319 

.2589 

.8232 

2.6178 

.5032 

.0074 

.3258 

1.0359 

3.2942 

Table  18 

Crude 

Estimates 

for  Network  II 

>=•5) 

II 

i 

.4931 

2 

.4973 

§ 

3 

.4947 

$ 

Exponential 

Service 

VAR 

WQ.S 

.0001 

.4800 

<.00005 

.3472 

Weibull  Service 


< .  00005 


Uniform 

Service 

< . 00005 

.2763 

.0002 

< .00005 

.0561 

< . 00005 

< .00005 

.0320 

< . 00005 

W 


WW 


Table  19 


Analytic  Estimates  for  Network  II  {Q 


Exponential  Service 

NODE 

RO.C 

VAR 

L 

RATIO 

1 

.4999 

< . 00005 

.  0005 

.0016 

2 

.4997 

<.00005 

.0064 

.0204 

3 

.4996 

< . 00005 

.0028 

.0090 

NODE 

WQ.C 

VAR 

L 

RATIO 

1 

.4829 

.0009 

.2191 

.6968 

2 

.3550 

.0004 

.1975 

.6280 

3 

.2299 

.0014 

.  7819 

2.4865 

Weibull  Service 

NODE 

RO.C 

VAR 

L 

RATIO 

1 

.4999 

<  .  00005 

.0012 

.0037 

2 

.4999 

<.00005 

.0049 

.0156 

3 

.4997 

<.00005 

.0028 

.0089 

NODE 

WQ.C 

VAR 

L 

RATIO 

1 

.3136 

.0005 

.3875 

1.2321 

2 

.1457 

.0001 

.2642 

.8402 

3 

.0995 

.0001 

.4986 

1.5854 

Uniform  Service 

NODE 

RO.C 

VAR 

L 

RATIO 

1 

.5000 

< .00005 

.0032 

.0101 

2 

.5000 

< .00005 

.  0040 

.0126 

3 

.  5000 

< .00005 

.0028 

.0089 

NODE 

WQ.C 

VAR 

L 

RATIO 

1 

.2759 

.0002 

.2534 

.8057 

2 

.0561 

< . 00005 

.3201 

1.0179 

3 

.0320 

<.00005 

.3655 

1.1624 

=  •5) 


U 

.0051 

.0649 

.0286 

U 

2.2158 

1.9970 

7.9071 


U 

.0118 

.0496 

.0283 

U 

3.9181 

2.6718 

5.0416 


U 

.0321 

.0401 

.0283 

U 

2.5621 

3.2369 

3.6964 


Table  20 


Modified  Estimates  for  Network  II  (9s- 5) 
Exponential  Service 


NODE 

RO  (M ) 

VAR 

L 

RATIO 

U 

1 

.4932 

.0001 

.2045 

.6503 

2.0680 

2 

.4981 

.0001 

.3856 

1.2263 

3.8996 

3 

.4950 

.0001 

.1539 

.4894 

1.5563 

NODE 

WQ(M) 

VAR 

L 

RATIO 

U 

1 

.4739 

.0009 

.2403 

.7643 

2.4305 

2 

.3413 

.0005 

.2904 

.9234 

2.9364 

3 

.2262 

.0007 

.3731 

1.1863 

3.7724 

Weibull  Service 


NODE 

RO  (M) 

VAR 

L 

RATIO 

U 

1 

.4957 

< .00005 

.1272 

.4045 

1.2863 

2 

.4976 

< .00005 

.2222 

.7065 

2.2467 

3 

.4961 

< . 00005 

.0743 

.2362 

.7511 

NODE 

WQ  (M  ) 

VAR 

L 

RATIO 

U 

1 

.3080 

.0003 

.2426 

.7716 

2.4537 

2 

.1448 

.0001 

.3964 

1.2604 

4.0081 

3 

.0971 

.0001 

.3494 

1.1110 

3.5330 

Uniform  Service 


NODE 

RO  (M) 

VAR 

L 

RATIO 

U 

1 

.4984 

< . 00005 

.0892 

.2837 

.9022 

2 

.4975 

< . 00005 

.0153 

.0486 

.1545 

3 

.4979 

<.00005 

.0205 

.0653 

.2077 

NODE 

WQ  (M ) 

VAR 

L 

RATIO 

U 

1 

.2473 

.0003 

.3616 

1.1498 

3.6564 

2 

.0554 

< .00005 

.  2763 

.8786 

2.7939 

3 

.0315 

< .00005 

.3875 

1.2324 

3.9190 

Table  21 


External  Estimates  for  Network  II  (P=.5) 
Weibull  Service 


NODE 

RO.C 

VAR 

L 

RATIO 

U 

1 

.  5004 

< . 00005 

.0236 

.0751 

.2388 

2 

.4990 

< . 00005 

.0792 

.2519 

.8010 

3 

.4996 

<.00005 

.0109 

.0346 

.1100 

NODE 

WQ.C 

VAR 

L 

RATIO 

U 

1 

.3225 

.0002 

.1824 

.5801 

1.8447 

2 

.1449 

< . 00005 

.0418 

.1330 

.4229 

3 

.1027 

<  . 00005 

.1038 

.3302 

1.0500 

Uniform  Service 


NODE 

RO.C 

VAR 

L 

RATIO 

U 

1 

.5002 

<.00005 

.3651 

1.1610 

3.6920 

2 

.4983 

< . 00005 

.  3060 

.9730 

3.0941 

3 

.4994 

<.00005 

.0821 

.2611 

.8303 

NODE 

WQ.C 

VAR 

L 

RATIO 

U 

1 

.2805 

.0004 

.5099 

1.6216 

5.1567 

2 

.0560 

<.00005 

.2915 

.9269 

2.9475 

3 

.0319 

< .00005 

.3700 

1.1766 

3.7416 

RESULTS  FOR  NETWORK  III 


Table  22  lists  the  steady  state  Jackson  values  for  the  uti¬ 
lization  factors  and  queue  times  of  Network  III  depicted  in 
Figure  7. 


Table  22 


Steady  State  Jackson 

Values  for 

Network  III 

B .  9 

Q- 

s.  5 

NODE 

RO  ( J ) 

WQ(J) 

RO  ( J ) 

WQ  ( J ) 

1 

.  9000 

6.0787 

.  5000 

.3750 

2 

.9001 

9.6031 

.  5000 

.4167 

3 

.9000 

15.1853 

.  5000 

.9377 

4 

.  9000 

5.7589 

.  5000 

.2500 

A  U  (O  H  ^WK)H  *.  Ul  IO  H 


Tables  23-26  list  the  results  for  Network  III  at  a  traf¬ 


fic  intensity  of  .9. 


Table  23 


Crude  Estimates  for  Network  III  ((>  =  .9) 


Exponential  Service 


NODE 


RO.S 

VAR 

WQ.S 

VAR 

.8912 

.0005 

6.0750 

7.5294 

•  8999 

.0003 

9.0400 

8.1951 

.8874 

.0006 

13.1571 

12.4821 

.8972 

.0003 

5.4977 

.5474 

Weibull 

Service 

.9320 

.0003 

6.5081 

6.2433 

.8943 

.0001 

4.3736 

.7746 

.8943 

.0003 

7.9743 

4.4438 

.8963 

.0001 

1.9755 

.0951 

Uniform 

Service 

.9353 

.0002 

6.0980 

7.2565 

.8939 

.0002 

2.8475 

.2466 

.8946 

.0004 

5.3534 

1.0206 

.8935 

.0002 

.5270 

.0056 

Tables  27-30  list  results  for  Network  III  at 


traffic 


intensity  of  . 5 


«5 

a 

f 


a 


u*< 
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|L  l 
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Table  24 


Analytic  Estimates  for  Network  III  (9«.9) 


Exponential  Service 


NODE 

RO.C 

VAR 

L 

RATIO 

U 

1 

.8989 

.0001 

.0382 

.1215 

.3864 

2 

.9053 

< .00005 

.0567 

.1804 

.5737 

3 

.9044 

.0001 

.0642 

.2041 

.6490 

4 

.9019 

.0001 

.1107 

.3519 

1.1190 

NODE 

WQ.C 

VAR 

L 

RATIO 

U 

1 

5.4466 

4.2202 

.1763 

.5605 

1.7824 

2 

7.8507 

4.6498 

.1784 

.5674 

1.8043 

3 

10.8947 

6.9876 

.1760 

.5598 

1.7802 

4 

5.0025 

.7746 

.4450 

1.4150 

4.4997 

Weibull  Service 


NODE 

RO.C 

VAR 

L 

RATIO 

U 

1 

.9132 

.0001 

.0606 

.1927 

.6128 

2 

.9007 

<.00005 

.0723 

.2300 

.7314 

3 

.9047 

.0001 

.0896 

.2850 

.9063 

4 

.8990 

<.00005 

.1021 

.3248 

1.0329 

NODE 

WQ.C 

VAR 

L 

RATIO 

U 

1 

5.7123 

4.1207 

.2075 

.6600 

2.0988 

2 

4.0249 

.5064 

.2056 

.6537 

2.0788 

3 

7.2919 

3.0851 

.2185 

.6948 

2.2095 

4 

1.8792 

.0878 

.2903 

.9233 

2.9361 

Uniform  Service 


NODE 

RO.C 

VAR 

L 

RATIO 

U 

1 

.9143 

.0001 

.1686 

.5361 

1.7048 

2 

.9019 

< .  00005 

.0261 

.0829 

.2636 

3 

.9043 

.0001 

.0464 

.1476 

.4694 

4 

.8977 

< . 00005 

.0772 

.2455 

.7807 

NODE 

WQ.C 

VAR 

L 

RATIO 

U 

1 

5.2239 

4.9088 

.2127 

.6765 

2.1513 

2 

2.5843 

.1150 

.1466 

.4663 

1.4828 

3 

4.5992 

.6964 

.2146 

.6823 

2.1697 

4 

.5125 

.0053 

.2977 

.9466 

3.0102 

& 
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Table  25 

Modified  Estimates  for  Network  III  ($g.9) 


Exponential  Service 


NODE 

RO  (M  ) 

VAR 

L 

RATIO 

U 

1 

.8927 

.0002 

.1575 

.5008 

1.5925 

2 

.9021 

.0001 

.1756 

.5583 

1.7754 

3 

.8927 

.0003 

.1860 

.5916 

1.8813 

4 

.8999 

.0003 

.2756 

.8763 

2.7866 

NODE 

WQ(M) 

VAR 

L 

RATIO 

U 

1 

5.7936 

8. 6035 

.3593 

1.1426 

3.6335 

2 

8.6201 

6. 8833 

.2641 

.8399 

2.6709 

3 

12.6169 

10.2152 

.2574 

.8184 

2.6029 

4 

5.3743 

.4023 

.2311 

.7350 

2.3373 

Weibull  Service 


NODE 

RO  (M ) 

VAR 

L 

RATIO 

U 

1 

.9340 

.0001 

.1216 

.3867 

1.2297 

2 

.8985 

<.00005 

.0726 

.2308 

.7339 

3 

.8973 

.0001 

.0997 

.3169 

1.0077 

4 

.8991 

.0001 

.1627 

.5173 

1.6450 

NODE 

WQ(M) 

VAR 

L 

RATIO 

U 

1 

6.1046 

5.8130 

.2926 

.9306 

2.9593 

2 

4.1846 

.6406 

.2601 

.8271 

2.6302 

3 

7.2838 

2.8174 

.1995 

.6345 

2.0177 

4 

1.9310 

.0734 

.2429 

.7724 

2.4562 

Uniform  Service 


NODE 

RO  (M  ) 

VAR 

L 

RATIO 

U 

1 

.9387 

.0001 

.0798 

.2539 

.8074 

2 

.8991 

< . 00005 

.0322 

.1024 

.  3256 

3 

.9015 

.0001 

.0516 

.1642 

.5222 

4 

.8961 

.0001 

.0945 

.3004 

.9553 

NODE 

WQ(M) 

VAR 

L 

RATIO 

U 

1 

5.6317 

4.3992 

.1907 

.6063 

1.9280 

2 

2.6042 

.1369 

.1746 

.5551 

1.7652 

3 

4.6425 

.6922 

.2133 

.6782 

2.1567 

4 

.5123 

.0054 

.  3068 

.9755 

3.1021 

Table  26 


External  Estimates  for  Network  III  (P-.9) 
Weibull  Service 


NODE 

RO.C 

VAR 

L 

RATIO 

U 

1 

.9374 

.0001 

.1164 

.3702 

1.1772 

2 

.8940 

.0001 

.2039 

.6483 

2.0616 

3 

.8977 

.0001 

.1371 

.4361 

1.3868 

4 

.8976 

.0001 

.1537 

.4888 

1.5544 

NODE 

WQ.C 

VAR 

L 

RATIO 

U 

1 

6.6310 

1.4742 

.0742 

.2361 

.7508 

2 

4.5533 

.7773 

.3156 

1.0036 

3.1914 

3 

8.2056 

4. 5019 

.3186 

1.0131 

3.2217 

4 

2.0323 

.0617 

.2041 

.6490 

2.0638 

Uniform  Service 


NODE 

RO.C 

VAR 

L 

RATIO 

U 

1 

.9393 

.0002 

.2631 

.8367 

2.6607 

2 

.8948 

.0002 

.2281 

.7255 

2.3071 

3 

.8969 

.0003 

.2488 

.7912 

2.5160 

4 

.8957 

.0002 

.3649 

1.1603 

3.6898 

NODE 

WQ.C 

VAR 

L 

RATIO 

U 

1 

6.0083 

2.5200 

.1177 

.3743 

1.1903 

2 

3.0071 

.3502 

.4466 

1.4202 

4.5162 

3 

5.4661 

1.0297 

.3173 

1.0089 

3.2083 

4 

.5298 

.0062 

.3503 

1.1141 

3.5428 

Table  27 


Crude  Estimates  for  Network  III  (p=.5) 
Exponential  Service 


NODE 

RO.S 

VAR 

WQ.S 

VAR 

1 

.4942 

.0002 

.3677 

.0009 

2 

.4969 

.0001 

.4379 

.0041 

3 

.4957 

.0001 

.9249 

.0074 

4 

.4976 

.0001 

.2400 

.0003 

Weibull  Service 

1 

.4970 

.0001 

.2373 

.0001 

2 

.4976 

.0001 

.2508 

.0007 

3 

.4983 

.0001 

.5540 

.0025 

4 

.4990 

.0001 

.1063 

<.00005 

Uniform  Service 

1 

.4976 

.0001 

.1867 

.0001 

2 

.4952 

.0001 

.1928 

.0003 

3 

.5007 

.0001 

.4486 

.0006 

4 

.4967 

.0001 

.0731 

<.00005 
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Table  28 

Analytic  Estimates  for  Network  III  (9=.5) 
Exponential  Service 


NODE 

1 

2 

3 

4 

NODE 

1 

2 

3 

4 


RO.C 

VAR 

L 

RATIO 

U 

4995 

<.00005 

.  0005 

.0015 

.0048 

5003 

< .  00005 

.0008 

.0025 

.0080 

5014 

<.00005 

.0019 

.0061 

.0194 

4996 

<  .  00005 

.0023 

.0073 

.0232 

WQ.C 

VAR 

L 

RATIO 

U 

3726 

.0003 

.0926 

.2945 

.9365 

4242 

.0027 

.2036 

.6476 

2.0594 

9284 

.0067 

.2839 

.9028 

2.8709 

2363 

.0001 

.1412 

.4490 

1.4278 

Weibull  Service 


NODE 

1 

2 

3 

4 

NODE 

1 

2 

3 

4 


RO.C 

VAR 

L 

RATIO 

U 

4995 

< .00005 

.0012 

.0038 

.0121 

5005 

<  .  00005 

.0007 

.0022 

.0070 

5013 

<.00005 

.0007 

.0022 

.0070 

4995 

<.00005 

.0012 

.0038 

.0121 

WQ.C 

VAR 

L 

RATIO 

U 

2390 

.0001 

.1336 

.4250 

1.3515 

2498 

.0005 

.2347 

.  7462 

2.3729 

5518 

.0011 

.1355 

.4310 

1.3706 

1055 

< . 00005 

.  1779 

.  5656 

1.7986 

Uniform  Service 


NODE 

1 

2 

3 

4 

NODE 

1 

2 

3 

4 


RO.C 

VAR 

L 

RATIO 

U 

4995 

<.00005 

.0007 

.0022 

.0070 

5002 

< . 00005 

.0005 

.0017 

.0054 

5015 

< . 00005 

.0005 

.0017 

.0054 

4996 

< .00005 

.  0002 

.0006 

.0019 

WQ.C 

VAR 

L 

RATIO 

U 

1872 

.0001 

.1245 

.3959 

1.2590 

1929 

.0001 

.  1413 

.4492 

1.4285 

4442 

.0003 

.1419 

.4513 

1.4351 

0733 

< . 00005 

.1222 

.  3886 

1.2357 

»-  J&t 
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Table  29 

V  y 

Modified 

Estimates 

for  Network  III 

(P  =  .5) 

M 

Exponential  Service 

jajy  NODE 

RO  (M ) 

VAR 

L 

RATIO 

U 

IP  1 

.4956 

.0001 

.0993 

.  3159 

1.0046 

vS  2 

.4992 

< . 00005 

.0154 

.0491 

.1561 

:£  3 

.4946 

< .00005 

.1109 

.3527 

1.1216 

SJj  4 

.4999 

.0001 

.2875 

.9142 

2.9072 

III  NODE 

WQ(M) 

VAR 

L 

RATIO 

U 

■1  1 

.  3673 

.0006 

.2097 

.6668 

2.1204 

&  2 

.4308 

.0023 

.1716 

.5458 

1.7356 

fe:  3 

.9149 

.0065 

.2779 

.8837 

2.8102 

.2413 

.0005 

.4626 

1.4712 

4.6784 

Weibull  Service 

!<J*  NODE 

m  i 

RO  (M ) 

VAR 

L 

RATIO 

U 

.4967 

< . 00005 

.0479 

.1522 

.4840 

>;•>'  2 

.4984 

< . 00005 

.0096 

.0306 

.0973 

5%J  3 

i 

K  NODE 

.4960 

< .00005 

.0430 

.1367 

.4347 

.4991 

< . 00005 

.0755 

.2402 

.7638 

WQ(M) 

VAR 

L 

RATIO 

U 

ft}  1 

.2354 

.0001 

.1958 

.6228 

1.9805 

.2479 

.0006 

.2623 

.8341 

2.6524 

3 

.5379 

.0013 

.1661 

.5283 

1.6800 

p  4 

.1052 

< .00005 

.3547 

1.1280 

3.5870 

Uniform  Service 

fr?  NODE 

RO  (M ) 

VAR 

L 

RATIO 

U 

1 

.4986 

< . 000G5 

.0039 

.0123 

.0391 

fi&|  2 

.4975 

< . 00005 

.0103 

.0326 

.1037 

m  3 

.4993 

< . 00005 

.0088 

.0281 

.0894 

pc  4 

.4974 

<  . 00005 

.0279 

.0886 

.2817 

f$V  NODE 

WQ(M) 

VAR 

L 

RATIO 

U 

E|  2 

.1862 

<.00005 

.1108 

.  3525 

1.1210 

.1893 

.0001 

.1394 

.4432 

1.4094 

.4400 

.0002 

.  1148 

.  3651 

1.1610 

I 

.0722 

< . 00005 

.0711 

.2260 

.7187 

1 

& 
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CHAPTER  V 


CONCLUSIONS 


The  purpose  of  this  research  was  to  study  the  application 
of  Jackson  networks  as  control  variates  in  queueing  simula¬ 
tions  in  order  to  make  some  general  conclusions  about  their 
effectiveness  for  variance  reduction.  These  conclusions 
hold  their  importance  in  that  they  add  to  the  store  of  pri¬ 
or  knowldge  an  analyst  can  draw  on  in  deciding  the  appro¬ 
priate  variance  reduction  technique.  Also,  these  results 
indicate  whether  or  not  continued  research  in  this  area  is 
warranted. 

The  results  of  this  study  indicate  the  potential  of  ana¬ 


lytic  controls  based  on  Jackson  networks  to  produce  vari¬ 
ance  reductions  in  utilization  factor  estimates.  Jackson 
based  controls  for  the  queue  time  estimates  were  not  as 


effective  as  the  utilization  controls  in  producing  variance 


reductions 


In  each  network  studied  the  queue  time  con¬ 


trol  variates  produced  little  or  no  variance  reduction  and 
indicated  the  potential  to  increase  this  estimate's  vari¬ 
ance.  In  some  cases  these  controls  could  increase  the  con¬ 


trol  estimate's  variance  up  to  eight  times  that  of  the 


crude  estimate's  variance.  The  analytic  controls  for  the 
utilization  factor  showed  more  promise. 

In  Network  I  the  analytic  controls  for  the  utilization 
factor  produced  variance  reductions  in  the  range  of  68  to 
88  percent  for  traffic  intensities  of  .9,  and  approximately 
99  percent  for  traffic  intensities  of  .5.  The  modified  ana¬ 
lytic  controls  produced  variance  reductions  of  approximate¬ 
ly  75  percent  only  for  the  uniform  service  time  case.  Per¬ 
formance  in  the  exponential  and  Weibull  cases  indicated  the 
potential  to  add  variance  to  the  estimate.  External  con¬ 
trols  for  the  utilization  factor  were  poor  and  again  indi¬ 
cated  the  potential  to  add  variance. 

Analytic  controls  for  the  utilization  factor  in  Network 
II  produced  consistent  variance  reductions  for  the  Weibull 
and  uniform  service  cases.  Reductions  for  these  controls 
ranged  from  80  to  94  percent  at  traffic  intensities  of  .9, 
and  approximately  98  percent  for  traffic  intensities  of  .5. 
The  modified  analytic  controls  performed  well  only  for  the 
uniform  service  case  at  traffic  intensities  of  .5.  Here  the 
reductions  ranged  from  71  to  93  percent.  External  controls 
were  generally  poor  with  the  exception  of  the  Weibull  ser¬ 
vice  case  at  the  .5  traffic  intensity  level  where  the  v  ri- 
ance  reductions  ranged  from  74  to  96  percent. 

Network  II  was  structured  so  that  each  service  node 
would  experience  the  same  effective  arrival  rate  in  steady 


75 


state  conditions.  The  service  rate  was  the  same  at  each 
node,  but  the  number  of  servers  was  varied  from  one  to 
three  (see  Figure  6).  The  purpose  was  to  observe  the 
impact  of  the  number  of  servers  on  control  variate  perform¬ 
ance.  The  results  did  not  indicate  any  observable  connec¬ 
tion  between  the  number  of  servers  and  control  variate  per¬ 
formance  . 

In  Network  III  the  analytic  controls  for  the  utilization 
factor  at  the  .9  traffic  intensity  level  showed  modest  per¬ 
formance.  Their  performance  at  the  .5  traffic  intensity 
level  was  greatly  improved  producing  variance  reductions 
of  approximatley  99  percent.  For  modified  analytic  con¬ 
trols  at  the  .9  traffic  intensity  level,  variance  reduc¬ 
tions  were  acceptable  only  in  the  uniform  service  case 
where  the  reductions  ranged  from  70  to  90  percent.  At  the 
.5  traffic  level  intensity  these  controls  showed  good  per¬ 
formance  for  both  the  Weibull  and  uniform  service  cases 
producing  reductions  in  the  range  of  76  to  98  percent. 
Performance  of  the  external  controls  was  generally  poor. 

The  results  did  show  the  analytic  control  variates  for 
the  utilization  factor  worked  well  at  the  .5  traffic  inten¬ 
sity  level.  The  same  statement  could  not  be  made  about  the 
queue  time  controls  since  their  performance  was  so  erratic 


at  both  traffic  intensity  levels  studied.  No  conclusive 
statement  could  be  made  concerning  the  impact  of  the  ser¬ 
vice  distribution  on  control  variate  performance. 
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The  Jackson  based  analytic  controls  indicate  promise  as 
effective  control  variates  for  the  utilization  factor.  One 
possible  explanation  for  the  difference  in  performance 
between  the  utilization  factor  and  queue  time  controls  may 
be  suggested  by  the  form  of  these  two  performance  measures. 
The  utilization  factor  is  a  ratio  based  on  effective  arri¬ 
val  and  service  rates.  The  queue  time  measure  has  a  more 
complex  form  incorporating  the  probability  distribution  of 
the  number  of  customers  at  a  service  node  and  the  fraction 
of  the  customer  load  carried  by  the  servers;  see  (8)  and 
(10).  This  may  suggest  the  variance  of  the  queue  time  or 
similar  measure  may  be  too  complex  to  be  fully  captured  by 
the  control  variate  approach.  This  should  not  preclude 
future  research  in  the  application  of  control  variates  in 
queueing  network  simulations.  One  approach  may  may  be  to 

^  A  A 

observe  different  forms  of  A  »  and  r.  to  obtain  the  con¬ 

trol  variates,  such  as  observing  p.(0)  of  (10)  directly 

1  \ 

from  the  simulation  rather  than  computing  it  from  A. »  JU.* 

A 

and  jl-  Another  approach  worth  considering  is  to  search  for 
models  which  provide  close  approximations  for  the  perform¬ 
ance  measures  of  interest  as  opposed  to  the  exact  analyt¬ 
ical  value  provided  by  the  Jackson  model.  The  use  of 
approximation  models  may  very  well  broaden  the  class  of 
queueing  networks  receptive  to  the  control  variate 
approach.  Whitt  [20]  has  investigated  the  use  of  open 


networks  to  approximate  the  performance  of  closed  queueing 
networks.  The  opportunity  of  expanding  the  use  of  control 
variates  in  open  systems  as  approximations  for  closed  net¬ 
works  would  be  enhanced  by  further  research  in  this  area. 

Another  possible  approach  to  improving  control  variate 
performance  is  to  obtain  a  more  precise  estimate  of  the 
control  coefficient  b.  One  way  to  accomplish  this  would  be 
to  increase  the  number  of  batches  in  a  replication.  A 
pilot  run  increasing  the  number  of  batches  from  25  to  50 
was  performed  on  Network  I  with  Weibull  service  at  the  .9 
traffic  intensity  level.  Little  improvement  was  noted  in 
the  performance  of  the  utilization  controls;  however,  con¬ 
siderable  improvement  was  seen  in  the  queue  time  controls. 
Variance  reductions  doubled  for  the  analytic,  modified  ana¬ 
lytic,  and  external  controls.  This  may  suggest  the  poor 
performance  of  the  Jackson  based  queue  time  controls  is  not 
solely  due  to  the  Jackson  model.  The  ability  or  inability 

it  , 

to  accurately  estimate  b  may  have  a  major  impact  on  con¬ 
trol  variate  performance  for  these  systems.  The  methodology 
for  estimating  the  control  coefficient  is  open  to  further 
study. 

Other  sources  which  may  explain  the  poor  performance  of 

the  control  variates  lie  in  the  methodology  of  this  study. 

*  * 

In  order  to  obtain  an  interval  estimate  and  a  value  for  b 


the  batch  means  approach  was  employed  rather  than  running  a 
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series  of  independent  replications.  This  is  not  an  uncommon 
practice  and  is  employed  to  reduce  the  cost  of  the  simula¬ 
tion.  The  batch  means  approach,  however,  produces  only 
approximate  independence  between  the  batches.  Determining 
the  batch  size  is  critical  to  this  independence  and  is  com¬ 
plicated  by  network  structure.  A  batch  size  may  work  well 
for  one  particular  node  and  not  as  well  for  the  remaining 
service  nodes.  Further  study  is  needed  to  determine  the 
usefulness  of  the  batch  means  approach  in  this  methodology. 
Another  source  for  error  is  the  initial  bias.  These  net¬ 
works  tended  to  have  long  and  erratic  initial  bias  periods. 
This  study  took  a  fairly  conservative  approach  in  deleting 
this  bias;  however,  further  study  of  the  initial  bias  in 
networks  is  needed  to  improve  the  application  of  control 
variates  for  steady  state  analyses. 

The  results  of  this  study  do  highlight  the  potential  of 
analytic  control  variates  in  simulation.  Depending  on  the 
parametric  model  selected  to  serve  as  the  basis  for  the 
control,  the  effort  required  to  obtain  variance  reduction 
would  be  6mall  compared  to  reducing  the  variance  through 
additional  run  time  or  the  second  simulation  required  by 
external  controls.  This  holds  considerable  promise  for 
automating  or  incorporating  the  analytic  control  approach 
in  existing  simulation  languages.  In  software  designed  for 
a  specific  user  this  approach  could  be  incorporated  by  the 


addition  of  a  statistical  collection  mechanism  and  a  rou¬ 
tine  to  derive  the  controls  from  these  statistics.  The  ben¬ 
efit  of  this  endeavor  would  be  to  avail  a  wider  range  of 
variance  reduction  techniques  to  the  user  community  and 
enhance  the  analysis  provided  through  computer  simulation. 


Appendix  A 
COMPUTER  CODE 


//  JOB  , 

✓/  RECION=768K, TIME=7 
/*JOBPARM  LINES=50O0,V=S,DISKIO=500O 
//  EXEC  SIM93CC,TIME.GO=6 
//CMP.SYSIN  DD  * 

’’ANTHONY  P.  SHARON  ADVISOR:  DR.  BARRY  L.  NELSON 
••DEPT:  ISE  THESIS  RESEARCH 

••APPLICATION  OF  JACKSON  NETWORKS  AS  EXTERNAL  CV 
• ’FOR  QUEUEING  SIMULATION  SYSTEM:  2  NODES 
• 'ARRIVAL: EXPONENTIAL  SERVICE: EXPONENTIAL 

’ ’BATCH  LENGTH:  1S0  INITIAL  DELETION: 10000 

’ ’NO.  OF  MACROS: 10  BATCHES  PER  MACRO:  25 

* ’ ATIMC I) 

’  ’BR(  I) 

*  I 

• ’BUSY! I) 

• ’CUST 

•  *  LAMBDA! I) 

' ’MU< II 

•  *  NODT 
”O.R(  I,J) 

• ’R( I.J) 

”RCUST(  I.J) 

”S(  I) 

• *STIM(  I) 

* ’ TCUST(  I) 

*  * 

’  ’ WTIM(  I) 

PREAMBLE  LAST  COLUMN  IS  72” 

EVENT  NOTICES  INCLUDE  RESET,  OUTPUT 
EVERY  ARRV  HAS  A  NODE. A 

DEFINE  RODE. A  AS  AN  INTECER  VARIABLE 
EVERY  EOS  HAS  A  CUST.E,  A  RODE . E 

DEFINE  CUST.E,  NODE . E  AS  INTECER  VARIABLES 
PERMANENT  ENTITIES 

EVERY  NODE  HAS  AN  ATIM,  A  BUSY,  A  LAMBDA,  A  MU,  AN  8, 

AN  AWQ,  AN  STIM,  A  TCUST,  A  VTIM  AND  OWNS  A  QUEUE 
DEFINE  BUSY,  S,  TCUST  AS  INTEGER  VARIABLES 
TEMPORARY  ENTITIES 

EVERY  CUST  HAS  AN  NODT  AND  MAY  BELONG  TO  THE  QUEUE 
DEFINE  RCUST  AS  A  2-DIMENSIONAL  INTEGER  ARRAY 
DEFINE  BR,  O.R,  R  AS  2-DIMENSIONAL  ARRAYS 
DEFINE  NM.C  AS  VARIABLES 
ACCUMULATE  A. BUSY  AS  THE  MEAN  OF  BUSY 
TALLY  A.WQ  AS  THE  MEAN  OF  WTIM 
TALLY  A.AR  AS  THE  MEAN  OF  ATIM 
TALLY  A. SR  AS  THE  MEAN  OF  STIM 
TALLY  C.AWQ  AS  THE  MEAN  OF  AWQ 
END  ’ ’ PREAMBLE 
MAIN 

DEFINE  I  AS  AN  INTEGER  VARIABLE 
LET  NM= 1 

CREATE  EVERY  NODE  (2)  ”NO.  OF  NODES 

RESERVE  BR(*,*> ,  R(*.*),  RCUSTC*,*),  O.R(*,*>  AS  2  BY  2 
READ  BUSY,  LAMBDA,  S 
START  NEW  RECORD 
READ  MU 

START  NEW  RECORD 
READ  BR 


INTERARRIVAL  TIME  AT  NODE  I 

INPUT  BRANCHING  PROBABILITIES  FROM 

NODE  I  TO  NODE  J 

NO.  OF  BUSY  SERVERS  AT  NODE  I 

CUSTOMER 

EXTERNAL  ARRIVAL  RATE  AT  NODE  I 

SERVICE  RATE  AT  NODE  I 

ENTRY  TIME  AT  A  NODE 

OBSERVED  BRANCHING  FROM  NODE  I  TO  J 

COMPARISON  ROUTING  MATRIX  BASED  ON  BR(I,J> 

NO.  OF  CUSTOMERS  ROUTED  FROM  NODE  I  TO  J 

NO.  OF  SERVERS  AT  RODE  I 

SERVICE  TIME  AT  NODE  I 

NO.  OF  CUSTOMERS  COMPLETING  SERVICE 

AT  NODE  I 

QUEUE  TIME  AT  RODE  I 
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FOR  1*1  TO  2  ,  DO  "HO.  OF  RODES 

FOR  J*  I  TO  2,  DO  "  RO . OF  RODES 
LET  CR*  BR( I ,  J)  ♦  CR 
LET  R( I,J)  =  CR 

LOOP 

LET  CR-e 

LOOP 

PR1RT  I  LIRE  THUS 

ECHO  IKPUT 

eriv  9  t lure 

FOR  1*1  TO  2,  DO  ’ ’RO.OF  RODES 

PRIRT  5  LIRES  WITH  I,  LAMBDA!  I),  HU(I>, 

BUSY!  I),  S! I)  THUS 
IRPUT  VALUES  FOR  BODE  * 

ARRIVAL  RATE:  *«.***« 

SERVICE  RATE:  **.**** 

RO.  BUSY  SERVERS:  ** 

RO.  OF  SERVERS:  ** 

SKIP  2  LIRES 

LOOP 

SKIP  2  LIRES 
LIST  BR 
6KIP  1  LIRE 
LIST  R 
SKIP  2  LIRES 

•• SCHEDULE  ARRIVAL  FOR  RODES  VI TH  EXTERNAL  ARRIVALS 
FOR  1*1  TO  2,  DO 
LET  RODE-  I 
LET  UA* RANDOM. F! NODE) 

LET  ATIM!  RODE)  = !  - 1 . 0 /LAMBDA!  RODE) ) *(  LOC.E.  F( UA)  ) 
SCHEDULE  AN  ARRV  GIVER  RODE  IR  ATIM! RODE)  UNITS 

LOOP 

SCHEDULE  A  RESET  1R  18800.8  UR ITS  "TIME  TO  DELETE  BIAS 
SCHEDULE  AN  OUTPUT  IR  10150.0  UNITS  "END  OF  FIRST  BATCH 
8TART  SIMULATION 
STOP 

END  "MAIN 

EVENT  RESET  "DELETES  BIAS,  RESETS  FOR  NEXT  BATCH 
FOR  EACH  NODE  RESET  THE  TOTALS  OF  ATIM.  STIM,  BUSY,  VTIH 
FOR  1*1  TO  2,  DO  "RO.  OF  MODES 

LET  ATIM!  I)*0 
LET  STIM!  I>*0 
LET  VTIM!  1)*8 
LET  TCUST!  I)  *8 

LOOP 

FOR  I* 1  TO  2,  DO 

FOR  J* I  TO  2, DO 

LET  RCUST! I,J>*0 
LET  O.R! I, J)=0 

LOOP 

LOOP 

RETURN 

END  "EVERT  RESET 

EVENT  ARRV  GIVEN  RODE  "EXTERNAL  ARRIVAL  AT  GIVER  RODE 

DEFINE  RODE  AS  AN  INTEGER  VARIABLE 

CREATE  A  CUST 

LET  ROOT! CUST)*  TIME.V 

LET  UA* RANDOM. F( RODE) 

LET  AT I M! RODE) * ( - 1 . 8/LAMBDA!  RODE) ) *( LOC . E . F( UA) ) 
SCHEDULE  AN  ARRV  GIVER  RODE  IR  ATIM!  RODE)  UNITS 
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IF  BUSY!  RODE) »  SI PODE),  PILE  COST  IP  QUEUE! PODE) 

ELSE 

LET  BUSY! PODE)*  BUSY! PODE)  *  1 

LET  WTIM!  PODE) *0 

LET  AWQ! RODE) *9 

LET  US* RANDOM. F! RODE+5) 

LET  STIM!  RODE)  * !  - 1 . •/'MU!  PODE)  ) *!  LOG.  E.  F(  US)  ) 

SCHEDULE  AP  EOS  GIVER  COST,  PODE  IP  STIM! PODE)  UP ITS 
ALWAYS 
RETURH 

ElfD  •  *  EVENT  AARV 

EVERT  EOS  GIVER  CUST,  PODE  "ERD  OF  SERVICE  EVERT 

DEFIRE  CUST,  RODE  AS  INTEGER  VARIABLES 

LET  TCUST!  RODE) *  TCUST(HODE)  ♦  1 

CALL  R0UTE2  GIVER  CUST, RODE 

RETURN 

ERD  *  ’  EVERT  EOS 
EVERT  OUTPUT 

DEFIRE  I.  J  AS  INTEGER  VARIABLES 
FOR  1*1  TO  2,  DO 

FOR  J>1  TO  2,  DO 

LET  0. R!  I ,  J) *  RCUST!  I,J)/TCUST(  I) 

LOOP 

WRITE  l.«/A.AR!I),  1.0/A.SRII),  0.R(I,1),  0.R(I,2) 

AS  4  D(ie,4)  US  IRC  UNIT  1 
WRITE  AS  /  USING  UNIT  1 

WRITE  A. BUSY! I)/S! I),  A.VQ(I)  AS  2  D! 10.4) 

USING  UNIT  1 
WRITE  AS  /  USING  UNIT  1 

LOOP 
LET  I*  I 
LET  C*C+I 

PRINT  6  LINES  WITH  RM.C,  I,  1+1, 1.0/A. AR!  I) ,  1.0/A. AR!  1+1)  , 

1 . 0/A. SR!  I )  ,  1 . 0/A.  SR!  1+ 1)  .  A.  BUSY!  D /S!  I)  ,  A. BUSY!  1+ 1  >/S(  D-l >  , 
A.WQII),  A.WQ!  I+I)  THUS 
RESULTS  FOR  MACRO  **.*  BATCH  **.* 

RODE  *  RODE  « 

ARRIVAL*  **.****  ARRIVAL*  **.**** 

SERVICE*  **.****  8ERVICE*  **.**** 

RHO*  *.*»**  RHO*  at.**** 

WQ=  ***.****  lfQ=  ***.**** 

PRINT  1  LIRE  THUS 

OBSERVED  ROUTING  PROBABILITY  MATRIX  !0.R) 

LIST  O.R 
SKIP  1  LIRE 
IF  C  LT  28.0 

SCHEDULE  A  RESET  ROW 
SCHEDULE  AP  OUTPUT  IP  180.0  UNITS 
ELSE  PRINT  1  LIRE  WITH  C.AWQ!1),  C.AWQ!2)  THUS 
OVERALL  HEAR  FOR  WQ!  I )***«.***«  WQC 2) =***.***» 

8KIP  1  LIRE 
LET  C*0 
LET  RM*RHM 

FOR  EACH  RODE  RESET  THE  TOTALS  OF  AWQ 
FOR  1*1  TO  2,  DO  "RO.  NODES 
LET  AWQ!  I)«0 

LOOP 

IF  RM  LE  10.0 

SCHEDULE  A  RESET  ROW 

SCHEDULE  AR  OUTPUT  IP  150.0  UNITS 


Y*viV^V,Ssy,,VJ**» 
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ELSE  STOP 
ALWAYS 
ALWAYS 
RE TURK 

END  ’  •  EVENT  OUTPUT 

ROUTINE  R0UTE2  GIVEN  CUST,  NODE  "FOR  TWO  NODE  NETWORK 
DEFINE  CUST,  NODE.  NEXT  AS  INTEGER  VARIABLES 
LET  DEST  *  RANDOM. Ft  8) 

IF  DEST  LE  R( NODE,  1 )  ,  LET  NEXT*  1 
LET  NODTt CUST) =TIME. V 

LET  RCUSTt  NODE , NEXT) *  RCUSTt NODE, NEXT)  +  i 
IF  BUSY! NEXT) *  S(NEXT),  FILE  CUST  IN  QUEUE ( NEXT) 

ELSE  LET  BUSYt  NEXT) *  BUSYC  NEXT)  *1 
LET  WTIMt  NEXT) =6 
LET  AWQCNEXT)** 

LET  US = RANDOM. Ft  NEXT+5) 

LET  STI M( NEXT)  =  ( - 1 . O/MUC NEXT) ) *( LOG. E. Ft  US) ) 

SCHEDULE  AN  EOS  GIVEN  CUST,  NEXT  IN  STI M( NEXT)  UNITS 

ALWAYS 

ELSE  IF  DEST  LE  R(N0DE,2) 

LET  NEXT  *  2 

LET  NODTt CUST)  «  TIME.V 

LET  RCUSTt  NODE, NEXT) *  RCUSTt NODE, NEXT)  +  1 
IF  BUSYt  NEXT)  *  S(NEXT) 

FILE  CUST  IN  QUEUEt NEXT) 

ELSE  LET  BUSYC  NEXT) =  BUSYt NEXT)  +  1 
LET  WTIMC  NEXT) *• 

LET  AWQt  NEXT) =• 

LET  US* RANDOM. Ft  NEXT+5) 

LET  STIMC NEXT)  =  t  - 1 . O/MUt  NEXT) ) *t  LOC. E. Ft  US) ) 

SCHEDULE  AN  EOS  CIVEN  CUST.  NEXT  IN  STIMC NEXT)  UNITS 

ALWAYS 

ELSE  DESTROY  CUST 
ALWAYS 

ALWAYS 

IF  QUEUEt NODE)  IS  EMPTY 

LET  BUSYt  NODE) -  BUSYt NODE)  -  1 

RETURN  _ 

ELSE  REMOVE  THE  FIRST  CUST  FROM  QUEUEt NODE) 

LET  WTIMt NODE)*  TIME.V  -  NODTt CUST) 

LET  AWQC  NODE) *  TIME. V- NODTt CUST) 

LET  US=RANDOM. Ft  NODE+5) 

LET  STIMC  NODE)  =  t - 1 . 0/MUC  NODE) ) *t  LOC. E. Ft  US) ) 

SCHEDULE  AN  EOS  CIVEN  CUST,  NODE  IN  STIMC NODE)  UNITS 

RETURN 

END  "ROUTINE  R0UTE2 

ROUTINE  SNAP. R 

LIST  TCUST 

SKIP  1  LINE 

LIST  RCUST 

SKIP  1  LINE 

LIST  ATTRIBUTES  OF  EACH  EOS  IN  EV. SCI. EOS) 

SKIP  1  LINE 

LIST  ATTRIBUTES  OF  EACH  ARRV  IN  EV.SCI.ARRV) 

RETURN 

END  "SNAP.R 
/* 


//  JOB  , 

//  REGI0N*768K 
/tJOBPARM  LINES* 5 00© 

//SI  EXEC  FORTVCC,  I I6L IB* SINGLE 
//FORT. SYS  IN  DD  * 

CCC  PROGRAM  ANALYZES  SIMULATION  RESULTS  FOR  VARIANCE  REDUCTION 
CCC  PROGRAM  COMPUTES  EFFECTIVE  ARRIVAL  RATES  AND  PERFORMANCE 
CCC  MEASURES  FOR  A  GIVEN  NETWORK  USING  AN  IKSL  ROUTINE 

CCC  TO  SOLVE  THE  BALANCE  EQUATIONS.  CRUDE  AND  CONTROL 

CCC  VARIATE  ESTIMATES  ARE  COMPUTED  AND  THEIR  VARIANCES 
CCC  ARE  COMPARED. 

CCC  DEFINITIONS  OF  MAJOR  VARIABLES 


CCC 

NN 

NO.  OF 

NODES 

ccc 

NB 

NO.  OF 

BATCHES  PER  MACRO 

ccc 

NM 

NO.  OF 

MACROS  PER  EXPERIMENT 

CCC  SUBSCRIPTS  FOR  DEFINITIONS 

CCC  UNO.  OF  NODES.  J«NO.  OF  BATCH,  K*NO.  OF  MACRO 
CCC  RO( I)  LONG  RUN  UTILIZATION,  NODE  I 

CCC  WQ(I)  LONG  RUN  QUEUE  TIME,  NODE  I 

CCC  NOTE:  STATISTICS  ARE  DEFINED  FOR  RO  ONLY;  NOTATION  IS 
CCC  SIMILAR  FOR  VQ 

CCC  ROJ( I)  JACKSON  STEADY  STATE  FOR  RO(I) 

CCC  ROA(l.J)  ANALYTIC  CONTROL  FOR  RO(I.J) 

CCC  ROS(I.J)  SIMULATION  ESTIMATE  FOR  RO(I,J> 

CCC  AROS(I.K)  MEAN  FOR  ROS(I,J>  OVER  MACRO  K 

CCC  AROA(I)  HEAR  FOR  ROA(I.J)  OVER  A  MACRO 

CCC  VROAC I)  VARIANCE  FOR  RO( I.J)  OVER  A  MACRO 

CCC  CROC  I )  COVARI ANCEC ROSC I.J). ROAC I.J))  OVER  A  MACRO 

CCC  BC ROC  I)  ESTIMATED  CONTROL  COEFFICIENT,  B,  FOR  ROC  I) 

CCC  BROCCI.K)  CONTROL  ESTIMATOR  FOR  ROC  I)  IN  MACRO  K 

CCC  ROCCI.K)  CONTOL  ESTIMATOR. B=  I ,  FOR  ROC  I)  IN  MACRO  K 

CCC  HROCCI)  MEAN  OF  ROCC I.K) 

CCC  VROCCI)  VARIANCE  OF  ROCC I.K) 

CCC  MB  ROCC  I)  MEAN  OF  BROCC  I ,  K) 

CCC  VBROCCI)  VARIANCE  OF  BROCCI.K) 

ccc  Haros ( n  mean  of  arosci.io 

CCC  VAROSCI)  VARIANCE  OF  AROSC I.K) 

CCC  VHROCCI)  VARIANCE  REDUCTION  FROM  ROCC I) 

CCC  VRBR0CC I)  VARIANCE  REDUCTION  FROM  BROCC I) 

CCC  MWQSC I.K)  OVERALL  MEAN  OF  WQS  AT  NODE  I.  MACRO  K 

CCC  CMWQSC I)  OVERALL  MEAN  OF  WQS  AT  NODE  I  FOR  ALL  K 

CCC  NOTE:  A  LIST  OF  OTHER  PROCRAM  VARIABLES  FOLLOWS 
CCC  AC  I, II)  EFFECTIVE  ARRIVAL  RATE  MATRIX 

CCC  BCI)  ARRIVALC  INPUT) /EFFECTIVE  ARRIVALC OUTPUT) 

CCC  IA  ROW  DIMENSION  OF  AC  I, II) 

CCC  1DCT  ACCURACY  TO  DECIMAL  PLACE  OF  LEQT1F  SOLUTION 

CCC  IER  LEQTIF  WARN  FLAC  C  ACCURACY  OR  SINGULARITY) 

CCC  LEQTIF  IMSL  LINEAR  EQUATION  SOLVER 

CCC  M  NO.  OF  RICHT  HAND  SIDES 

CCC  N  NO.  OF  ROWS  IN  B(  I) 

CCC  R(i.Ii)  PR.  OF  ROUTING  FROM  NODE  I  TO  II 

CCC  S(I)  NO.  OF  SERVERS  AT  NODE  I 

CCC  WKAREA(I)  DIMENSION  CT  OR  EQ  TO  N 

CCC  MU( I)  SERVICE  RATE  AT  NODE  I 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
CCC  DECLARE  AND  DIMENSION  VARIABLES 

INTEGER  NN,  NM,  NB,  M,  N,  IER,  IA.  IDCT.  S  .FACT, CRUDE 
REAL  HWQS(2, 1«) ,CMWQS(2) ,CVWQS(2) 

REAL  RO J ( 2 ) ,  R0A(2,50) ,  AROS(2,10),  AROA(2),  VR0AI2). 


♦  CROC 2) ,  BCROC 2) ,  BROC(2,10),  ROC! 2, 10),  MROC! 2) , 

♦  VROCC2),  MB ROC ( 2) ,  VBROCC2),  MAROSI2),  VAR0SC2), 

♦  VRROCC2) .  VRBROCC2) 

REAL  WQJI2),  WOAI2.50),  AWQSI2.10),  AWQAC 2) ,  VWQAC  2) , 

+  CWOI2),  BCWQC2),  BVQCC  2, 10) .  WQC(2.10>,  MWQCC2), 

+  VWQC(2>,  MBWQCC  2) ,  VBWQCC2),  VRWQCC2),  VRBWQCC2) 

REAL  B(2>.  WKAREAC 4) ,  MU(2),  RC2.2),  A(2,2) 

CCC  LABELED  COMMON  STITT 

COMMON  /SIMMEA/  ROSI2.25) ,  VQSC2.25),  8(2) 

CCC  READ  IN  NETWORK  PARAMETERS  AND  IMSL  ARGUMENTS 
R£AD(5,*)NN,  NM.  NB.  M.  N,  IA.  1DGT 
CCC  READ  JACKSON  PARAMETERS  AND  FORM  A(I.II)  MATRIX 
DO  20  Is 1 , NN 

READ! 5,*)Bi  I)  ,  MU(  I)  <  S(I),  R(I,I),  R(I,2) 

20  CONTINUE 

CCC  READ  OVERALL  MACRO  MEANS  FOR  WO 
DO  2S  K= 1 , NM 

READ(5,*)MWOS( 1,K) ,  MWQSC2.K) 

23  CONTINUE 

DO  40  I = 1 , NN 

DO  30  1 1  *  1 , NN 

IF(  I . EQ. 1 1) THEN 

A(  1 ,  1 1 )  *  I-R(  1 ,  1 1 ) 

viev 

A(  1 , 1 1 )  *  -R<  1 1 , 1 ) 

END  IF 
30  CONTINUE 

40  CONTINUE 
CCC  ECHO  INPUT 

WRITE( 6  41) 

41  FORMAT! 5 0',I3X,  ’ECHO  INPUT’) 

WRITE! 6 , *) ’ NN* ’ «  NN,  ’NB=  * ,  NB,  ’HM=  * ,  NM.  ’H=’,  M,  *N=’.  N, 

+  * IA= ’ ,  IA,  ’ IDGT= ’ .  I DOT 

DO  60  I* 1 .NN 

WRITE! 6,42) I ,  B!I),  MU!  I),  S!  I) 

42  FORMAT! ’0’, ’NODE’,  I3.2X, ’B! I ) = * .F10.4.2X, 'MU! I) * ’ ,FI0.4, 

♦  2X, ’SI  1)  =  ’ , 13) 

DO  80  1 1 r 1 , NN 

WRITE! 6, 43) I,  II,  RII.II),  A! I, II) 

43  FORMAT! ’0* , ’PROBABILITY  FROM’ , 13, IX, ’TO’ , 13, IX, ’«’ ,F10. 

+  IX,  ’A(  I,  ID*  ’  ,F10.4> 

50  CONTINUE 
60  CONTINUE 

CCC  SOLVE  FOR  JACKSON  EFFECTIVE  ARRIVAL  RATES 
CCC  CALL  IMSL  ROUTINE  LEOT1F 

CALL  LEGT1F! A,  M,  N,  IA,  B,  IDCT,  WKAREA,  IER) 

IF! IER.CT.0)THEN 

WRITE! 6,*) ’ERROR  FLAG= ’ ,  IER,  ’JACKSON’ 

STOP 
END  IF 

CCC  COMPUTE  JACKSON  MEASURES 
J«  1 

CALL  PERF1NN,  NB.  J,  B.  MU,  ROA,  WUA,  CRUDE) 

DO  70  I* 1 , NN 

ROJ! I) -ROA! I.J) 

WQJ!  D-WQA!  I.J) 

70  CONTINUE 

CCC  ECHO  JACKSON  MEASURES 
DO  110  I-I.NN 

WRITEC6.75) I,  ROJ!  I),  WQJ! I) 
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75  FORMAT< *0’ , ’ JACHSON  VALDES  NODE’ , 13, 2X,  *R0= ’ ,Fie.4,2X, 

+  'WQ=’,F10.4> 

lie  CONTINUE 

CCC  FOR  EACH  MACRO  COMPUTE  MEASURES  AMD  CVS 
DO  270  K= 1 , NM 

CCC  FOR  EACH  BATCH  READ  READ  PARAMETERS  AND  SIM.  ESTIMATES 
DO  140  J-l.NB 

CCC  FOR  EACH  NODE  READ  PARAMETERS  AND  ESTIMATES 
DO  130  I* 1 , NN 

READ(3,*)B(  I)  ,  MU(I>,  R(I,1),  R(I,2> 
READ<3.*)R0S( I, J) ,  WOS(I.J) 

130  CONTINUE 

DO  136  1*1, NN 

DO  135  1 1 = 1 , NN 

I  Ft  I.Ett.  IDTHEN 

A(  I, II)*1-R(  I, II) 

vi  or 

A(  1, 1 1  >  =  -R( II, I) 

END  IF 

135  CONTINUE 

136  CONTINUE 

CCC  COMPUTE  EFFECTIVE  ARRIVAL  RATE  FOR  A  BATCH 
CCC  CALL  1MSL  ROUTINE  LEttTIF 

CALL  LEQT1F( A,  M,  H,  IA,  B,  IDCT,  WKAREA,  IER) 

IF(  IER. CT.0) THEN 

WRITE(  6,*) ’ERROR  FLAG  BATCH’ ,  J,K,  IER 

STOP 

END  IF 

CCC  COMPUTE  BATCH  MEASURES  FOR  EACH  NODE 

CALL  PERFtNN,  NB,  J.  B,  MU,  ROA,  VQA,  CRUDE) 

140  CONTINUE 

NRITEC6, 14 1)K.  CRUDE 

141  FORMAT! *0’ , ' IN  MACRO’ ,  I3.2X, ’RHO  GE  1 .0’ , I3.2X, ’TIMES’ ) 
CCC  COMPUTE  MEANS,  VAR’S,  COV,  BSTAR,  FOR  BATCH  OUTPUT 

DO  240  1=1. NN 

DO  150  J*1,NB 

AROS(  I ,  IO  =  ROS( I , J)+AROS(  I.K) 

AWQS( I , K) =  WQS(  I,J)+AWQS( I.K) 

AROA(  I ) =  ROA(  I, J)+AROA<  I) 

VROACI)*  ROA(  I, J)*ROA( I,J)+VROA( I) 

AWQA(I)=  WQA( I , J)+AWQA( I) 

VWOAC  I )  =  VQA(  I, J)*VttA(  I,  J)+WQA(  I) 

150  CONTINUE 

CCC  COMPUTE  MEAN,  VAR  FOR  MEASURE  AT  NODE  I 

AROS(  I , K) =  AROS(  I ,K)/(NB) 

VROA(  I)  *VROA(  I )  /'(  NB-  1 ) -AROA(  I)*AROA(  I)-' 

+  ( <  NB- 1 ) *( NB) ) 

AROA(I)*  AROA(I)/(NB) 

AVQS(  I , K) =  AVOS!  I , K) /( NB) 

WQA(I)*  VWOA(  I)/(NB-1)-AVQA(  I)*AWQA(  I)/ 

*  1 1  NB- 1 )*(  NB) ) 

AWQA( I ) =  AWQA(I)/(NB) 

CCC  COMPUTE  COVARIANCES 

DO  160  J* I , NB 

CRO( I)*  (ROS( I, J)-AROS( I,K))«(ROA( I,J)-AROA( I))+ 
+  CRO(I) 

CWtt( I)=  (VQ S(  I,J)-AVQS( I , K) ) *( VQA( I,J)-AKQA(  I))  + 

*  CVQ(  I ) 

160  CONTINUE 

CCC  COMPUTE  COV  AND  BSTAR 
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CRO!  I)«CR0!  I)/!HB-1> 

CVQC  I>*  CVQC  I)/(NB-1> 

}  BCROC  !)■  CROC I>/VROA( I) 

BCVQC  I ) ■  CVQC IJ/WQA! I) 

CCC  COMPUTE  CONTROL  VARIATE  ESTIMATES  FOR  RUN  K 

BROC(  1 ,K) *  AROS( I,IO -BCROC I)*(AROA( I)-ROJ! I>) 

ROC( I,K)*AROSC I,K)-AROA( I)+ROJ( I) 

BWQCt  I,E)«  MWQS(  I,K>-BCVa(  I)*!AVQA!  I)-VQJ!  D) 

VQCC I ,K)«  HVQS! I,K)-AVQAC I)+WOJ( I) 

CCC  PRINT  DATA  SUMMARY  FOR  RUN  K 

VR1TE!6,170)I,K 

170  FORMAT ( ’O’ ,  ’RESULTS  FOR  NODE* .  I3.5X,  ’MACRO’,  13) 

VRITE!6.180)R0J! I).  AROSCl.K),  CRO(I>.  BCROC  I)  .BROCC  I,X>  , 
♦  ROCCI.K) 

180  FORMAT! '  *,  ’ROJ*’,  F10.4.2X, • ARCS* • ,  F10.4.2X, 

+  ’COV(RO)*’ ,FI0.4,2X, ’ BSTAR= • .F10.4.2X, ’BROC* ’ .FI0.4, 

4  2X, *ROC* ’ ,  FIS. 4) 

WIUTEC  0,200)  WOJC  I)  ,  HVU3C  I,K)  ,  CWU  I)  ,  BCWUC  I)  .BVQCC I, K)  , 
4  VQCC I ,  K) 

200  FORMAT!  •  *VQJ*  *  .F1S.4.2X.  •  AVQS*  *  .F10.4.2X,  'COV!UQ>*  • , 

4  F1S.4.2X,  ’BSTAR*  '  .F10.4.2X,  ’BVOC*  *  .F10.4.2X,  ’WOO  »  ,F1S.4> 

240  CONTINUE 

CCCC  INITIALIZE  BATCH  MEASURE  ARRAYS 
DO  260  I* I ,NN 
DO  2S0  J*1,NB 
ROA! I.J)*0 
VQA! I ,  J)*0 

280  CONTINUE 

AROA!  I>«0 
VROA!  I) *0 
CRO!I)*0 
BCROC I)*0 
AVQAC 1)«0 
VWQA!  I>*0 
CVQ! I)>0 
BCVQC  I)*0 

260  CONTINUE 

CRUDE* 0 
270  CONTINUE 

CCC  COMPUTE  VARIANCE  REDUCTION  SUMS 
DO  390  I-l.NN 

DO  280  K*1,NH 

HR0C!I>*  ROC! 1 , K) +KROC!  I) 

VROCiD*  BOCC I,K)*ROCC  I.D4VR0C!  I) 

MVQCC I>*  VQCC 1,K)4MVQCC I) 

WQC!D*  VQCC  1,K)*VQCC  I,K>4WQCC  I) 

MBROC!  D*  BROCC 1 , K>  +MBROCC  I ) 

VBROCC I) *  BROCC I ,K)*BR0C! I , Kl+VBROC! I) 

MB VOCC I)*  BVQCC 1 ,K)4  MB VOC! I ) 

VBVQCC I)»  BVQCC I,K)*BVQC! 1,K>4  VBVQC! I) 

HAROS C I)*  AR0S!1,K)4  HAROS! 1) 

VAROS!!)*  AROS! I,K)*AROS! I,K)4  VAROS! I) 

CHVQS!  I)*  HWQS!I,K)4  CMVQS!  I) 

CWQSi  I) *  HNQS!I,K)*HMQS!I,K>4  CVVQSC  I) 

280  CONTINUE 

CCCCC  COMPUTE  VARIANCE  /MEANS 

VROCC  I)*  VROCC I )/C NM-D-MROCC  I)*MROC!  1 ) /! ! NH- 1 ) *NH) 

MROCC 1)>  MROCC I)/NM 

WQCC  I>*  VVOCC  I ) /! NM“  1 ) -MVQCC  I)*HVQC!  I)/t !NN-1)*NH) 

MVQC!  I>*  MVQCC  D/NM 
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ccc 


ccc 


290 

300 


310 


320 


330 


340 


VBROC!  I )  ■  VBROC ( D/INH- l)-HBROC!  I)*HBROC!  I)/<!KM-1)*NM) 
MBROC! 1 ) «  MBROC! I)/NM 

VBWQC!  I) -VBWQC!  I)/!  HM-D-MBVQC!  I)*HBVGC!  I)X!1NH-1)*NM) 
MBVQC! !)■  HBWQC!I)/NM 

VAROS(  D-VAROS!  D/TNM-D-MAROS!  I)*MAROS!  1)S!!NM-1)*NH) 
HAROS!  I)*  KAROS! II/NM 

CVWQS!  I)»CVVGS1  I)/!NH-I)-GMVQS!  I)*CHVQS!  1 ) ✓( ( TO- 1 ) *NH) 
CMWQS!I>-  CMWQS!  I) /TIM 
COMPUTE  VARIANCE  REDUCTION 

VRROC! I)-! VAR06! I)-VROC( IO/VAROS! I) 

VRBROC( I)-!VAROS! I)-VBROC( I) >/VAROSl I) 

VRWQC!  I)>(GWQS(  II-WQC!  1))/CVVQS!  1) 

VRBWQC(  D-ICVWQS!  I) -VBWQC!  I) > /CVWQS!  I) 

PRINT  VARIANCE  RDUCTIONS 
WRITE! 6,290) 1 

FORMAT! *0* •  'VARIANCE  REDUCTION  AT  NODE  ’ ,  13) 

WRITE! 6 ,300) ROJ!  I)  ,  HAROS!  1),  VAROS!  I) 

FORMAT!'  '.  'ROJ»  ’ ,  F10.4.3X,  'MEAN  SIM  RO*  * ,  F10.4.3X, 
'VAR! SIM  RO)-',  FI0.4) 

WRITE!4,310)MBR0C! I),  VBROC! I),  HROC! I) , VROC! I) 

FORMAT!'  'MEAN  BROC- ' ,  F10.4.3X,  ’ VAR! BROC) « • ,  F10.4, 

3X,  'MEAN  ROC-',  F10.4.3X,  'VAR!ROC>-',  F10.4) 

WRITE! 6 , 320) VRBROC! I > ,  VRROC!I) 

FORMAT! '0'.  'VAR  REDUCE  BROC-',  F10. 4. 5X, * VAR  REDUCE  ROC* ' , 
F10.4) 

WRITE! 6. 330) WQJ! I).  CMWQS!I>,  CVWQS! I) 

FORMAT! '0' ,  'WQJ-',  F10.4.3X,  'MEAN  SIM  WO**,  F10.4.3X, 

’ VAR! SIM  WO)-*,  F10.4) 

WRITE! 6 , 340) MB WOC ( I ) ,  VBWQC!!),  RWOC!I),  VWQC! I) 


FORMAT! 


MEAN  BWOC- ' ,  F10.4.3X, ' VAR!BWQC)* ’ ,  F10.4, 


330 

390 


3X, 'MEAN  WOC-’,  F10.4.3X, 
WRITE! 6,330) VRBWOC! I) ,  VRWOC! I) 
FORMAT! 


VAR!  WOC)  ■ 


F10.4) 


CCC 


800 


CCC 


CCC 


CCC 

CCC 


.  'VAR  REDUCE  BWOC-',  F10.4,8X. 'VAR  REDUCEWOC-' 
►  F10.4) 

CONTINUE 

STOP 

END 

INTEGER  FUNCTION  TO  COMPUTE  A  FACTORIAL 
INTEGER  FUNCTION  FACT! ISERV) 

ESERV- I 

IF!  ISERV. GT.  DIMEN 
DO  800  I- 1, ISERV 
ISERV- KSERV* 1 
CONTINUE 
END  IF 

FACT* KSERV 

RETURN 

END 

SUBROUTINE  TO  COMPUTE  PERFORMANCE  MEASURES 
SUBROUTINE  PERF!NN,  NB,  J,  B,  HU,  ROA,  WQA, CRUDE) 

COMMON  /S IMMEA/ROS! 2,28),  WQS!2,25),  8(2) 

INTEGER  NN.NB.J,  FACT,  8,  CRUDE 
REAL  B! 2) ,  MU!2>,  ROA! 2, 25),  WOA!2,28),  LOA 
COMPUTE  MEASURES  FOR  EACH  NODE 
DO  820  I- I, NR 
TI-0 

ROA!  I,J)-B(  I)/!S!  I)*MUi  D) 

IF  RO  GREATER  THAN  ONE,  REPLACE  RO  WITH  .9999 
TO  COMPUTE  MEASURES 

IF! ROA!  I, J)  .CE.  OTHER 


V 

«' 

I* 

4‘ 


<!t 

& 

A 

»*’ 


,S 

$ 

to 


R0A<  I,J)».W» 

CRUDE*  CRUDE* I 

SiTwm  l.n*K  i»«s<  l»/<MCT>»<  '»« 

SIR  CORTIRUE 

WQA(I,J>*  UtA/Bd) 

820  CORTIRUE 
RETURN 
CRD 

J/coffrwreSi  Sd  wii.'is»M.'raBi.»i»»*<®1 
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